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A COMPUTER PROGRAM TO CALCULATE 
RADIATING VISCOUS STAGNATION STREAMLINE 
FLOW WITH STRONG BLOWING 


By G. Louis Smith and L. Bernard Garrett 
Langley Research Center 

SUMMARY 

A computer program (program LEE) has been developed to calculate the fully 
coupled solution of the radiating viscous stagnation streamline flow with strong blowing. 
Program LEE was developed for the study reported in NASA TR R-388. The present 
report describes the digital computer program, including FORTRAN IV listing, flow 
charts, instructions for the user, and a test case with input and output. Program LEE 
is available through COSMIC. 

INTRODUCTION 

An implicit finite difference solution to the viscous shock-layer stagnation stream- 
line flow with radiation and strong blowing is presented in references 1 and 2. Program 
LEE (Langley computer program A3899) was developed to make the calculations reported 
in references 1 and 2. A detailed radiative transport code (RATRAP) that accounts for 
the important radiative exchange processes for gaseous mixtures is used (ref. 3). 
Chemical equilibrium of the flow is assumed, and the equilibrium composition of the air- 
injection products mixture is computed by a free-energy minimization routine (FEMP) 
within RATRAP. 

The program may be used for the stagnation streamline of a two-dimensional body, 
an axisymmetric body, or a three-dimensional body with two orthogonal planes of sym- 
metry. The difference between these cases is mainly in the continuity equation. For the 
three-dimensional case, there are two components of velocity around the body to consider, 
but the momentum equation is the same for both components. 

A description of the computer program is presented in this paper. Included are 
program listing, flow charts, instructions for the user, and a test case with input and 
output listings. Copies of the program LEE computer deck are available through COSMIC 
(Computer Software Management and Information Center at Barrow Hall, University of 
Georgia, Athens, Georgia). 



SYMBOLS 


a 


b 


Cp 


H 


h 


I 

K 


k 

N 


N 


Pr 


N 


Re 


PL 


‘R,y 




tangential velocity gradient, 

ox 

tangential velocity gradient, 

oz 

specific heat at constant pressure 
binary diffusion coefficient 
specific total enthalpy 
specific static enthalpy 
index for nodal point 
body curvature 

body curvature in z -direction for nonaxisymmetric body 
thermal conductivity 

number of node points in finite difference solution 
Prandtl number 
Reynolds number 
static pressure 

2 -5 

free-stream static pressure, dynes/cm (where 1 dyne = 10” newton) 
net radiative heat flux in y -direction 
radius of body, cm 
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radius measured from axis of symmetry of body (see fig. 1) 
free -stream velocity, cm/sec 

component of velocity parallel to surface in x-direction 
component of velocity normal to surface 
component of velocity parallel to surface in z -direction 
distance measured along body surface (see fig. 1) 
distance measured normal to body surface (see fig. 1) 

u . i. ^ 

distance measured along body surface in plane normal to plane in which x 
is measured 

mass fraction of ith chemical species 
mass fraction of foreign or ablation products 


x=0 


z=0 


-'n 


transformed shock -layer thickness, / p (y) dy 

'O 

error criterion for density 
transformed coordinate normal to body surface (see fig. 2) 


't/ 

•/n 


P (y) dy 


scale factor, 1 + Ky 



K. 

1 


dummy scale factor 


K scale factor for curvature in z-direction, 1 + K^V 

z 

M viscosity 

p density 

Subscripts: 

C carbon 

H hydrogen 

N nitrogen 

O oxygen 

n nth nodal point; n = 1 at the body ( i? = 0 ) and n = N at the shock (17=1) 

s shock 

00 free -stream conditions 

Superscript: 

1 ith iteration 

Unprimed symbols are nondimensional. All dimensional symbols are primed. 
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Figure 2. - Flow -field coordinate system in the transformed coordinates. 
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PROBLEM DESCRIPTION 

The equations which govern radiating viscous flow along the stagnation streamline 
of a blunt body in nondimensional form are as follows (ref. 1 or 2); 


Continuity: 


— (kpv) = - Kpa 


d ( K ^pv) _ 


dy 


= -2 k pa 


d(KKzPV) 


dy 


x-momentum: 




(Two dimensional) (la) 


(Axisymmetric) (lb) 


(Three dimensional) (Ic) 


(2a) 


z -momentum: 


N. 


Re 


_^/ db\ 
dy ( dy j 


+ K^pbv= a ^ 


y -momentum: 


dv 


dy 


(2b) 


(3) 


Energy: 


1 ^ /^i'/^ dH\ 

''Re \^Pr 


- \KK.PV-^ 
^ dy 


1 d / i dv\ d / 
^Re dy(^Np/dyj^dy 


KK . Q 



(4) 
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Diffusion: 


d 

dy 


<cKiPDi2 



K/CPV 


d Op 

"d^ 


= 0 


(5) 


These equations are written in the coordinate system shown in figure 1. In these equations, 
K = 1 for the two-dimensional case, = k for the axisymmetric case, and for 

the three-dimensional case. 


For the case of a three-dimensional body with two orthogonal planes of symmetry 

X is the distance from the stagnation point around the body in one plane of symmetry and 

z is the distance measured in the other plane of symmetry. There are two corresponding 

orthogonal components of velocity u and w parallel to the body surface. The stagnation 

streamline equations, which are the limit forms of the fluid-flow equations, involve the 

terms a =-^^and b = For two-dimensional flow, b and z -curvature vanish. For 
3 x az 

axisymmetric flow, b is identical with a. The quantities have been nondimensionalized 
as follows (where primes denote dimensional quantities): 


x,y,r 


x’,y',r' 

^b 


u,v,w 


u’,v’, w* 
U’ 




K = K’ 




h,H = 


h'.H' 


(uL)' 


Di2 = 


N 


p’ u’ r: 

oo oo - O 


Re " uL 


P- = 


P' 

P'„ 


y p' At' \3 






N 


c' u' 

p ^ 


Pr - k- 


P = 


P’ 


a = 


b = 




k’ 


k' 

s 



3u’ 

U’ 

OO 

t 

3x 

% 

9w' 


u' 

^ oo 


( 6 ) 


The symbol Rj|^ which is used for nondimensionalization is not restricted to being the nose 
radius of curvature. For example, for a flat disk, R^ can be the disk radius. The 
following restrictions are made on the problem: 
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(1) The gas is in local thermodynamic and chemical equilibrium. 

(2) Diffusion is governed by a binary diffusion model. 

(3) Radiative energy transport occurs within a one -dimensional, infinite, planar 
slab (tangent-slab approximation). 

(4) The gas density must be high enough that the shock wave is thin relative to 
the shock layer. 


The transform variable 


V 



p(y) dy 


(7) 


is defined where 



P(y) dy 


( 8 ) 


It is seen that this transformation is defined such that the shock wave is fixed at ?? = 1. 
The transformation of equation (7) is applied to equations (1) to (5) so that they can be 
written, respectively, as follows: 


Continuity: 


d (Kpv) 
d T) 


5«a 


(Two dimensional) (9a) 


d (k^pv) 
d T} 


2 Sks. 


(Axisymmetric) (9b) 


d (KKz pv) = j^a = 3/c b 

dTj ^ 


(Three dimensional) (9c) 


X- momentum: 

KP A/pp ^ 

dr}[ dvj s dv 

^^Re 


- pa + Kp av = -/? 


(10a) 
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z -momentum: 


/ X 2 

d / db\ ^ 



+ K pbv = - 13 
z z 


y -momentum: 


drj 



Energy: 


(10b) 


( 11 ) 


. lOC.Pfl 

0 J, 

dr, Np^ 




KK. pV 
1 


df) 


"ANpr 


PV 



+ 5N^ 4- 

Re 


d_ 

dr) 




(12) 


Diffusion: 


drjl' 


' 2„ '*“f' 

““i" °12 dT, 


+ 


Skk^PV 



(13) 


The user is cautioned that the normal velocity v has been redefined as negative in the 
positive y- or r, -direction in the aforementioned equations. The purpose of this sign 
change is to make v positive across the shock layer, except within the layer of injected 
products. The flow -field coordinate system is shown in figure 2. The boundary condition 
at the wall for the diffusion equation has been generalized from that used in references 1 
and 2 for strong blowing and is 


’’12 


(1= 0) 


This boundary condition is appropriate also for moderate blowing rates. 
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Equations (10), (12), and (13) are written infinite difference form by using a 
central -difference formula for second derivatives and windward differencing for first 
derivatives in convective terms. Equations (9) and (11) are integrated for (kk^pv) and p, 
respectively, subject to boundary conditions at the shock wave. Up tb 100 node points 
may be used for the solution. These node points are evenly spaced ini? . 

PROGRAM ORGANIZATION 

The governing equations (9) to (13) are solved by an iterative scheme as shown in 
figure 3 and are discussed in reference 2. This section of the report describes subpro- 
grams used to carry out the computation. A flow chart of each subprogram is included 
except where computation is straightforward. A full listing of program LEE is given in 
appendix A. The input is described in appendix B, and the output is described in appen- 
dix C. Langley library subroutines which are used are described in appendixes D, E, 
and F. A list of each subprogram and its description is given as follows: 


Subprogram 

LEE 

GETREDY 

INITIAL 

RANHUG 

FOFXRO 

UPDATE 

FLOW 

CHEXRO 

MASS 

Y MOMNTM 
SPECIES 


Description 

Initializes program, iterates density and 
flow equations, and prints results 

Reads spectral data input for RATRAP 

Initializes flow variables 

Solves Rankine -Hugoniot relations for post- 
shock conditions 

Function used by RANHUG 

Option of input for initial profiles 

Iterates flow equations to solution 

Computes density, compares with previous 
value 

Solves continuity equation for v 

Solves y -momentum equation for p 

Solves diffusion equation for and elemental 
mass fractions 


10 


























Subprogram 


Description 


X MOMNTM 

ENERGY 

RADARAY 

REFTEMP 

D2DP2 

DIMPROP 

RHOSY 

DIFCOEF 

VISPRAN 

VISCOS 

PRANDTL 

DIFT 

TINT 

TRIDIAG 


Solves x-momentum equation for a and, if 
necessary, for b 

Solves energy equation for h and H 

Calls RATRAP for radiation heat flux at NIC 
points and linearly interpolates for remain- 
der of N points; also updates density 

Computes post -shock temperature 
2 

Computes ]3 = - ^ P - 
ax^ 

Computes dimensional p' and h', given 
nondimensional p and h 

Equation of state, p (p, h, a 


Computes binary diffusion coefficient 

Computes Prandtl number and nondimensional 
viscosity profiles, given p and h profiles 

Computes dimensional viscosity using tables 
of reference 4 

Computes Prandtl number using tables of 
reference 4 

Given a profile y^^Cx^, computes profile 


Given a profile y^(x^, computes profile 


y( O d t , where J is a dummy variable 


Solves Ax = D by Potter's method where A 
is a tridiagonal matrix and x and D are 
vectors 
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Subprogram 


Description 


SET Sets one array equal to another 

ITEST Computes number of C(I) differing from 

B(I), where C and B are arrays, by 
amount exceeding a specified E 

RATRAP Radiation transport code (ref. 3) 

The radiation transport code RATRAP, developed by K. H. Wilson (ref. 3), is 
used as a subroutine. It uses a number of subprograms which are listed as part of 
appendix A; these are all described in reference 3. Included among these is a free-energy 
minimization program (FEMP) which solves for the mass fractions of the individual 
species. The many spectral data required for RATRAP are read in by GETREDY at the 
beginning of the computations. Radiative flux is computed at each of a desired number N 
of evenly spaced points. The N is input in NAMELIST NAM4. 

Langley library subroutines of the SCOPE 3. 0 operating system which are used 
in the LEE program are ITRl, FTLUP, and DISCOT. These are described in appen- 
dixes D, E, and F, respectively. 

In some of the flow charts, the operations "damp and clip” appear. "Damp" refers 
to computing a value for a given variable by weighting the value most recently computed 
with its value from the previous iteration, for each nodal point. For example, the value 
computed for density would be 


P 


n 




(i-1) 

n 


+ (1 - dp) P 


(i) 

n 


where d^ is the damping factor for density. The purpose of damping is to suppress 
numerical oscillations from one iteration to the next, and thereby improve convergence. 
The selection of dp is discussed in references 1 and 2. "Clip" refers to the require- 
ment that a profile not go below or above some value. For example, an intermediate 
iteration may give an excessive radiation divergence term, thus causing an unrealistic 
dip in the enthalpy profile. This condition is alleviated by requiring that H(I) be at least 
as large as the wall value. When this adjustment is required, a message is printed. 
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PRIMARY SUBPROGRAMS 
Program LEE 

LEE is the main program. It calls subroutines GETREDY and INITIAL, which 
input data and initialize profiles, respectively. Then, LEE directs the solution of the 
fluid -flow equations and the density computation by successive iteration until convergence 
is achieved. The flow chart for LEE is given as follows: 
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Subroutine INITIAL 


The INITIAL subroutine performs those tasks which must be done before the 
interaction of the profiles can begin. The conditions of the problem are read in. Post- 
shock conditions are computed. Profiles of various flow quantities are initialized by one 
of two options. If lUPDATE = 0, linear and second-order profiles are assumed initially 
across the shock layer. If lUPDATE > 0, the UPDATE subroutine is called to input 
profiles. The flow chart for INITIAL is given as follows: 
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Subroutine RANHUG 


The RANHUG subroutine computes the properties behind a normal shock wave. 

An initial guess of pressure and static enthalpy is made and an initial guess for density 
ROSHOCK is computed on this basis. Starting with this value, a Newton -Raphson sub- 
routine ITRl iterates to compute ROSHOCK to the desired accuracy. The ITRl subroutine 
is described in appendix D. The flow chart for RANHUG is given as follows: 











Function Subprogram FOFXRO 


The FOFXRO function subprogram calls ITRl to compute ROSHOCK. It uses 
the Rankine -Hugoniot equations and equation of state (as represented by RHOSY) to 
evaluate density ROSHOCK. The flow chart for FOFXRO is given as follows; 








Subroutine UPDATE 

The UPDATE subroutine reads NAMELISTs for the profiles required to begin 
the iteration for the solution, provided that lUPDATE >0, This option allows the user 
to begin the iteration with the best possible guess. The flow chart for UPDATE is given 
as follows: 
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Subroutine FLOW 


The FLOW subroutine iterates the fluid-flow equations by computing the a(r?), 
v(t?), h{r)), p(ij), and Op(r?) profiles and shock -wave standoff distance for a given p(t?) 
pr;ofile. The flow chart for FLOW is given as follows: 
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Subroutine CHEXRO 


The CHEXRO subroutine computes a density profile RO(I) from the most recent 
p(I), h(I), and elemental profiles. This RO(I) profile is compared with the previous 
profile and if the accuracy criterion is satisfied, ICHXRO is set equal to 0; otherwise, 
ICHXRO is set equal to 1. The RO(I) profile is then damped by using the previous value. 
Subroutine ICHXRO is used by LEE to determine whether the solution has converged or 
should be continued. The flow chart for CHEXRO is given as follows: 
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Subroutine MASS 


The MASS subroutine integrates the tangential flow gradient pa(T? ) to compute 
the normal mass flow pv (rj) in accordance with the continuity equation. Also, the shock- 
layer thickness and the transformation from n to y are computed. Finally, the physical 
thickness of the shock layer DELBA is computed. There are three geometry options; 
IDIM = 1 for two-dimensional flow, IDIM = 2 for axisymmetric flow, and IDIM = 3 for 
the stagnation streamline of a three-dimensional body. 

When the computation enters the MASS subroutine, the geometry option flag IDIM 
is checked, in effect selecting which of equations 9(a), (b), or (c) is to be used. The 
right-hand side of the selected equation is computed and integrated, without the factor 
of 5. The appropriate scale factors are then divided out, and S is computed so that 
(pv) = 1. The remaining computations follow in sequence. The flow chart for MASS 
is given as follows: 
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Subroutine Y MOMNTM 


The Y MOMNTM subroutine integrates the y -momentum equation by using a 
density p(I) and velocity profile v(I). The flow chart for Y MOMNTM is given as follows: 
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Subroutine SPECIES 


The SPECIES subroutine solves the binary diffusion equation for the mass fraction 

profile a j,(I) of the injected mass. The elemental mass fractions a ^^(I) are then computed 

from a_(I), where i = C, H, N, and O. The flow chart for SPECIES is given as follows: 

F 
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Subroutine X MOMNTM 


The X MOMNTM subroutine solves the viscous x-momentum equation for the 
tangential velocity gradient profile a(I). Because this equation is nonlinear in a(I), a 
quasilinear approach is used, whereby the solution is computed by iteration. The flow 
chart for X MOMNTM is given as follows: 
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Subroutine ENERGY 


The ENERGY subroutine solves the energy equation. The RATRAP code is used 

9Qr ^ 

to compute radiation heat fluxes ^(1), from which — ^ is computed. The energy 
equation is solved for H(I), after which h(I) is computed. The flow chart for ENERGY 
is given as follows; 
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Subroutine RADARAY 


The RADARAY subroutine calls RATRAP for the radiation heat fluxes QR(I) at 
NIC points and linearly interpolates for the values of QR(I) (where QR(I) = q„ (I)) at the 

R5 y 

remaining points if IRFLUX >0. If ERFLUX g 0, QR(I) is set equal to zero. Also, density 
RO(I) is updated according to the most recent p(I), h(I), and elemental profiles. If damping 
of the density profile RO(I) is desired, RODAMP is set equal to some positive integer in 
the NAMELIST NAM4 input. If so, the old density profile, stored in XNRO(I), is used. 

The flow chart for RADARAY is given as follows: 
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GAS -PROPERTIES SUBPROGRAMS 


Several of the gas -properties subprograms are for specific gases, for example, 
air. For other gases, these subprograms should be replaced. The VISCOS and PRANDTL 
function subprograms select values by using tables from reference 4, which constitute 
most of the BLOCK DATA subprogram. The DIFCOEF subroutine computes a binary 
coefficient for a mixture of atomic hydrogen and atomic nitrogen. 

The RATRAP subroutine is used to compute the equilibrium composition, density, 
and temperature of a gaseous mixture of hydrogen, oxygen, nitrogen, and carbon by use 
of a free -energy minimization subroutine (FEMP). 

. USAGE 

Program Information and Deck Configuration 

Program LEE was written in the FORTRAN IV language for the Control Data 6600 
computer system under the SCOPE 3. 0 operating system. A field length of 150 000_ 

O 

locations is required. Computing time depends on the blowing rate and the accuracy of 
the initial profiles assumed; but it is of the order of 20 to 30 minutes for cases with radi- 
ation for 21 node points and radiation computed at 11 node points. The vast majority of 
the tirne is spent in the very complicated radiation computations. The flow -computation 
time for fixed density and no radiation is approximately 1 second for each iteration with- 
in the FLOW subroutine, so that 2 to 4 seconds are required for a converged solution. 

The computation time using FEMP for computing the equilibrium chemistry of the air- 
ablation product mixture and without radiation computations was found to be 1^, 2^, 
and 3^ minutes for nondimensional blowing rates of 0, -0.1, and -0.2, respectively. 

In this case, most of the time is spent in chemistry computations. The following sketch 
shows the deck configuration needed for execution: 
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Input Description 

The RATRAP data deck is supplied with a permanent code and is not normally 
modified. A listing of the input cards is included with the sample case. Following the 
RATRAP data deck is the LEE case -input data, which is loaded by NAMELIST. The 
first group of data, in NAM4, are numbers which may be varied at will but are usually 
left unchanged. The next group of data, in NAM5, specifies the case of interest. The 
final data are optional; these are the profiles which may be read in, arid they are in 
NAMl, NAM2, and NAM3. The input symbols are nondimensional except where indicated. 
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The NAM4 data 
N 

IDIM dimension option: 

IRODAMP 

HDAMP 

RODAMP 

EP 

EROV 

EA 

EPSX 

ERO 

MAXTIME 

NIC 

CURV 

CURVZ 

BS 

D2PDZ 

The NAM5 data 
IRFLUX 


are given as follows: 

number of points in flow-field solution, up to 100 

(1) Plane flow 

(2) Axisymmetric flow 

(3) Three -dimensional flow 

option for damping density: 0 no damping used on density 

1 damping used on density 

damping used on static enthalpy profile 

damping used on density profile 

accuracy criterion for pressure 

accuracy criterion for pv 

accuracy criterion for tangential velocity gradient, a ; 

d'X 

also used for b = — 

dZ 

accuracy criterion used in solving x-momentum equation 

relative accuracy criterion for density 

maximum number of iterations allowed within FLOW 

number of points for which RATRAP computes radiative fluxes, 
up to 100 

curvature of body at x = 0 in plane in which x is measured 

curvature of body for three-dimensional case at x = 0, z = 0 
in plane in which z is measured (used only for IDIM = 3) 

(used only for IDIM = 3) 

(used only for IDIM == 3) 
are given as follows: 

radiation option: 0 no radiation computed 

1 radiation computed 


dz /x=0, z=0 
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lUPDATE 


ALPINJC 

ALPINJO 

ALPINJH 

ALPINJN 


option for inputing initial guess profiles: 

0 no profile read in, assume linear 

1 read in profile 

mass fraction of carbon in ablation products injected 
mass fraction of oxygen in ablation products injected 
mass fraction of hydrogen in ablation products injected 
mass fraction of nitrogen in ablation products injected 


2 


D2DP2C 

^ A- coefficient 

AS 

value of a at shock wave 

PINP 

2 

free -stream pressure, dynes/cm (where 1 dyne = 10 

RHOINP 

free -stream density, g/cm^ 

UINP 

free -stream velocity, cm/sec 

RB 

body radius at x = 0, cm 

ROVW 

mass injection rate (negative) 

HW 

wall static enthalpy 

ALPINFN 

free -stream mass fraction of nitrogen 

ALPINFO 

■ y 

free -stream mass fraction of oxygen 

ALPINFC 

free -stream mass fraction of carbon 

ALPINFH 

free -stream mass fraction of hydrogen 


The NAMl data are given as follows: 

RO(I) density profile, initial guess 

ROV(I) profile of mass flow normal to body, initial guess 

The NAM2 data are given as follows: 

A(I) profile of tangential velocity gradient, initial guess 

H(I) static enthalpy profile, initial guess 


newton) 
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P(I) 

DELTA 


pressure profile, initial guess 

initial guess of transformed shock -layer thickness 


The NAM3 data are given as follows: 

AFOR(I) profile of mass fraction of injected species where 

I = 1, 2, ... N 

A sample set of input data is given in appendix B. This case corresponds to a 
304.8-cm (10-ft) nose radius body at an altitude of 60.96 km (200 000 ft) with a velocity 
of 15. 25 km/sec (50 000 ft/sec). This is t37pical of a manned reentry for a Mars return. 
Air is being injected at a high rate (0. 1 of free -stream value of mass flux). 

Discussion of Output 

The output of program LEE consists of printing only. Output for the sample case 
is shown in appendix C. The first page of output, printed by INITIAL, gives the conditions 
of the case, followed by the postshock conditions. The NAMELIST inputs NAM4 and 
NAM5 are printed on following pages, also by INITIAL. Next, equilibrium composition 
profiles as computed by FEMP from the initial pressure and static enthalpy profiles are 
written by RATRAP. If the radiation option is used, radiative heat fluxes are written. 

The next page shows the initial flow profiles, printed by INITIAL. For each pass through 
FLOW, the profiles of the various flow parameters are written. These are followed by 
profiles of elemental distribution, printed by SPECIES. If the radiation option is used, 
radiative heat fluxes are printed by 'RATRAP. Once the iteration for the full -coupled 
solution has converged, the message "VARIABLE DENSITY SOLN CONVERGED" is 
printed by LEE. The output symbols are nondimensional, except where indicated. 

The symbols for postshock conditions are given as follows: 

PSHOCK static pressure 

ROSHOCK density 

VSHOCK velocity 

HSHOCK static enthalpy 
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The symbols for flow profiles are given as follows: 


ETA 

n 

Y 

y 

QR 

radiative heat flux 

CAPH 

total enthalpy 

V 

V 

EMU 

M 

PRAN 

Ntj 

Pr 

YOYS 


ROV 

pv 

A 

a 

P 

p, pressure 

RO 

p, density (latest value) 

XNRO 

p, density (previous value) 

H 

h, static enthalpy (latest value) 

XNH 

h, static enthalpy (previous value) 

T 

T', temperature, kelvin 

AMEC 

a^, elemental mass fraction of carbon 

AMEN 

a elemental mass fraction of nitrogen 

AMEO 

a Q, elemental mass fraction of oxygen 

AMEH 

a elemental mass fraction of hydrogen 

AEOR 

a total mass fraction of injected products 

r 


AONC 


total mass fraction of free -stream gas 



The symbols for radiative output are given as follows: 


TOTAL QMINUS 
TOTAL QPLUS 
QRAD 


total radiative heat flux toward shock wave due to atomic 
2 

emission, W/cm 

total radiative heat flux toward body due to atomic line 
2 

emission, W/cm 

2 

net radiative heat flux, W/cm 


Error Messages 

The following table gives the error diagnostics as the message and the corre- 
sponding subprogram which prints the message: 


Message 

Diagnosis 

Subprogram 

FLOW SOLUTION NOT CONVERGED 

Fatal 

FLOW 

MOMENTUM FAILS TO CONVERGE 

Fatal 

X MOMNTM 

FEMP BLEW 

Fatal 

RATRAP 

MAXIMUM ITERATIONS EXCEEDED IN RANHUG .... 

Fatal 

RANHUG 

ZERO DERIVATIVE IN RANHUG 

Fatal 

RANHUG 

H(I) LESS THAN H(l) 

Nonfatal 

ENERGY 

AFORI(I) LESS THAN l.E-06 

Nonfatal 

SPECIES 

AFOR(I) GREATER THAN 0.999999 

Nonfatal 

SPECIES 

TEMPERATURE TOO LARGE, T' = ? 

Nonfatal 

VISCOS 

TEMPERATURE TOO LARGE, T’ - ? 

Nonfatal 

PRANDTL 

T DID NOT CONVERGE, T' - ? 

Nonfatal 

FEMP 

N-R DID NOT CONVERGE 

Nonfatal 

FEMP 


If the flow solution does not converge, the enthalpy and density profiles should be 
checked for a sequence of iterations. Adjustment of RODAMP and HDAMP may be indi- 
cated. Also, an improved approximation for the initial profiles reduces the number of 
iterations required. The momentum solution is very reliable. If it should fail to con- 
verge, an improved initial profile should be an adequate fix. 

Experience with FEMP has shown that the message FEMP BLEW is usually due 
to negative mass fractions or imrealistically low enthalpies. These problems are 
alleviated in the computations by the CLIP operation. 
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No difficulties have been experienced with the computation of postshock conditions. 
However, diagnostic messages are included in the RANHUG subroutine to provide for 
this contingency. For some conditions, the Newton -Raph son (N-R) procedure in FEMP 
has difficulty converging. In this case, the message "N-R DID NOT CONVERGE" is 
printed and the computation is reinitiated. 

The nonfatal error messages are to alert the user to items of which he should be 
aware, but which do not disrupt the computation. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., August 7, 1973. 
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APPENDIX A 


PROGRAM LISTING 

PROGRAM LEE( INPUT, OUTPUT, TAPE5=fNPUT, rAPE6=0UTPUTI A I 

A 2 

LAVER WITH A 3 

equilibrium chemistry and a a 

EQUILIBRIUM RAOIATICN A 5 

A 6 

COMMON /CA/ AllOn, 8(1011 A 7 

COMMON /CY/ Y( 1011 ,Y0YSI 1011 A 8 

COMMON /CRO/ ROIlOll A 9 

COMMON /CETA/ ETAIlOll A 10 

COMMON /CV/ VIlOll A 11 

COMMON /CP/ PIlOll A 12 

COMMON /CROV/ ROVllOll A 13 

COMMON /CEMU/ EMUIlOll A 14 

COMMON /CCAPH/ CAPHIlOll A 15 

COMMON /CH/ HIlOll A 16 

COMMON /CPRAN/ PRANIlOll A 17 

COMMON /CXRO/ XNROIlOll A 18 

COMMON /CXH/ XNHIlOll A 19 

COMMON /CQR/ QRdOll A 20 

COMMON /CT/ TdOll A 21 

COMMON /CAMFI/ AMEC dOO 1 , AMEN ( 100 1 , AMFO dOO 1 , AMFHd 00 1 A 22 

COMMON /CFORONC/ AFOR ( 20 1 1 , AONC I 20 1 1 A 23 

COMMON /MISCL/ N, IDIM,CURV,DETA, DELTA, PSHOCK,REY,HH,EP,EROV,EA, MAX A 24 

1TIME,MLDIFPT,NIC,ER0 A 25 

COMMON /DAMP/ IR00AMP,R00AMP,H0AMP A 26 

COMMON /RADCCMP/ DELBA,IRFLUX A 27 

COMMON /REFVALL/ EMUSHOK A 28 

CALL GETREDY A 29 

CALL INITIAL A 30 

CALL FLOW A 31 

CALL CHEXRO dCHXROl A 32 

IF (ICHXRC.EQ.il GO TO 1 A 33 

WRITE (6,51 A 34 

WRITE (6,61 DELTA, OELBA,REY A 35 

WRITE (6,71 A 36 

DO 2 1=1, N A 37 

WRITE (6,81 ETA( I I ,V( I 1,QR(I 1,CAPH( Il,Vd l,EMU( I l,PRAN( I l,YOYSd 1 A 38 

CONTINUE A 39 

WRITE (6,91 A 40 

DO 3 1=1, N A 41 

WRITE (6,101 ETA( I l,RQV( I I ,A( I l,P( II,ROd 1,XNR0( I 1,H( I I, XNH( II A 42 

CONTINUE A 43 

WRITE (6,111 A 44 

DO 4 1=1, N A 45 

WRITE (6,101 ETA(ll,T(II,AMFC(Il,AMFN(II,AMFO(Il,AMFH(II,AFORdl,A A 46 

lONCdl A 47 

CONTINUE A 48 

CALL PRINTIT A 49 

STOP A 50 

A 51 
A 52 

FORMAT (1H0,5X,31HVAR1 ABLE DENSITY SCLN CONVERGEDI A 53 

FORMAT dHl,5X,6HDELTA=E16.8,2X,6H0ELBA=E16.8,2X,4HREY»E16.81 A 54 

FORMAT (4X ,6HETA( 1 l,9X,4HY( I 1 ,10X,5hCRl I 1,8X,7HCAPH( I 1, 10X,4HV(I 1 , A 55 

19X,6HEMU( 1 1 ,8X,7HPRAN( I 1 ,8X,7HYOVS(I I I A 56 

8 FORMAT (8E14.6I A 57 

9 FORMAT (4X,6HETAdl,8X,6HROV(II,9X,4HA(Il,10X,4HP(Il,10X,5HROdl,8 A 58 

1X,7HXNR0( I I ,9X,4HH( I 1 ,9X,6HXNH(I 1 1 A 59 

10 FORMAT (8E14.61 A 60 

11 FORMAT (4X,6HETA( I 1 ,9X,4HTd I ,9X,7HAMFCd 1 ,8X,7HAMFN( I I ,8X,7HAMF0( A 61 

II I ,6X,7HAMFH( I I ,7X,7HAF0R( II ,7X,7HACNC( 1 1 1 A 62 

END A 63- 
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APPENDIX A — Continued 


1 

2 


5 


SUBROUTINE INITIAL B 

CQNMCN /CA/ A(IOX) ,B(101 I B 

CGf'MON /CBETA/ BETAIlOl* B 

COMMON /CTEMP/ TEMP ( 1 01 ) , T IMP (101 » B 

COMMON /CY/ Y(101),Y0YStl01) 6 

COMMON /CETA/ ETA(lOl) B 

COMMON /CP/ PdOll B 

COMMON /CRC/ BO(lOl) B 

COMMON /CROV/ ROV(lOl) B 

COMMON /CSCH/ SCHIlOlIfSCKlOlI 8 

COMMON /CV/ V(101» 8 

COMMON /CH/ H(lOl) B 

COMMON /CT/ T(lOl) B 

COMMON /CFORONC/ AFOR ( 20 U , AONC ( 20 U B 

COMMON /FRSTRM/ H IN , P 1 N , RHOI NP , U I NP , R 6 B 

COMMON /CAMFI/ AMFC ( 100 ) , AMFN ( 100 » , AMFO ( 100 1 1 AHFHI 100 ) B 

COMMON /MFINF/ ALP INFC t ALP INFN , ALP INFO , ALP INFH B 

COMMON /MFINJ/ ALP I N JC t A LP I N JN, ALP I N J 0 , ALP IN JH 8 

COMMON /RANK/ ROSHOCK .VSHOCK , HSHOCK , T SHOCK B 

COMMON /MISCL/ Nf IOIM,CURV,DETA, DELTA, PSHOCK,REY,HW,EP,EROV*EA, MAX 6 
1TIME,MLD1FPT,N1C,ER0 b 

COMMON /REFVALU/ EMUSHOK B 

COMMON /RADCQMP/ DELBA,IRFLUX B 

COMMON /DAMP/ IROOAMP.RQOAMP.HDAMP B 

COMMON /CD2DP2/ 020P2C 8 

COMMON /C3D/ CURVZ,BS,BETAZ(101I B 

NAMELIST /NAMAZ N , I 01 M , I ROOAMP , HDAMP , ROOAMP ,EP , EROV, EA , EPSX, ERO , HA B 
1XTIME,NIC,CURV,CURVZ,BS,02P0Z B 

NAMELIST /NAM5/ I RFLUX , I UPDATE , ALP I NJC, ALPINJO, ALP IN JH, ALP IN JN,020 B 
1P2C,AS,PINP,RH0INP,UINP,RB,R0VW,HW,ALPINFN,ALPINF0,ALPINFC,ALPINFH B 
READ (5,NAMA) $REA0( 5 ,NAM5 I B 

WRITE (6,20t B 

IF ( IRFLUX.GT.OI WRITE (6,191 B 

WRITE (6,171 R8,UINP,RHOINP,AS,D2DP2C,PINP,ROVW,ALPINJH,ALPINJC, AL B 
IPINJN, ALPINJO B 

WRITE (6,181 ALPINFN,ALPINFO,ALPINFC,ALPINFH B 

PIN=PINP/(RHOINP*UINP*UINPI B 

HIN=3.5*PIN*0.5 B 

CALL RANHUG B 

WRITE (6,NAMAI $WRITE ( 6 , NAM5 I B 

PEY=RHOINP*UINP*RB/EMUSHCK B 

RO(NI=RCSHOCK B 

XN=N-1 B 

CETA=1./XN B 

R0V(N)=1. B 

A(1)=0. B 

TEMP=AS-A(1I B 

ETA(11=C. B 

DO 1 1=2, N B 

ETAd ) = ETA( I-1)*0ETA B 

DO 2 1=1. N 8 

BETA! n=02DP2(ETA( I) 1 B 

IF ( I0IM.NE.31 GO TO A B 

DO 3 1=1, N B 

BETAZd l=02PDZ B 

B(I) = BS*ETA(n B 

CONTINUE 8 

IF ( lUPDATE.GT.0 1 GO TO lA B 

0ELTA=1. B 

0ELBA=.05*RB B 

DO 5 1=1, N B 

ROV(II=(ROV(N)-ROV(1)I*ETA(I)*ROV(1) B 

ROV( I l=ROV(l l*(l.-ROV( 11 1*ETA( I 1*ETA(II B 

H( I 1=(HSH0CK-HW1*ETA( 1 1*HW B 

Am = A(ldTEMP*ETA(Il B 

CONTINUE B 

00 6 1=2, N B 


1 

2 

3 

A 

5 

6 

7 

8 
9 

10 

11 

12 

13 

lA 

15 

16 

17 

18 

19 

20 
21 
22 
23 
2A 

25 

26 

27 

28 

29 

30 

31 

32 

33 
3A 

35 

36 

37 

38 

39 
AO 
A1 
A2 
A3 
AA 
A5 
A6 
A7 
A8 
A9 

50 

51 

52 

53 
5A 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 
67 
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APPENDIX A — Continued 


6 ROU )=RCSHOCK/ETA(H B 68 

RQ(1'»=RQ(2I*5. B 69 

DO 7 I=lfN B 70 

FGREIGN=0. B 71 

IF IROV(I).GT..O) FOREIGN=l. B 72 

AFOR(I)=l. -FOREIGN B 73 

7 AONCm=FOREIGN B 74 

CO 8 I=lf N B 75 

AMFCm = AFORI I )*ALPINJC + ACNCm*ALP !NFC B 76 

AMFNm=AFORm*ALPINJN+AONCm*ALPINFN _ B 77 

AMFOI I )=AFOR( n*ALPINJOFAONC( n*ALPINFO B 78 

AMFH(H=AFOR( I)*ALPINJH+AONCm*ALPINFH B 79 

8 CONTINUE B 80 

CO 9 1=1, N B 81 

VI I l=ROV( I l/RC< n B 82 

SCHm = l.+DELTA*ETA( n/R0( I I B 83 

SCim = l.+DELTA*ETA(n/RC( I) B 84 

9 CONTINUE B 85 

TEPP=-.5*V(NI B 86 

PSTAG=PSHOCK-TEPP B 87 

00 10 1=1, N B 88 

10 P( n=PSTAG+TEMP=ETAlII*=2 B 89 

DO 11 1=1, N B 90 

11 TEMPI H=1./R0( n B 91 

CALL TINT <N,OETA,TEMP,T.IMP» B 92 

DO 12 1=1, N B 93 

vm=DELTA=TIMP( I ) B 94 

SCHI I )=1.+CURV=0ELTA*T IMP! n 6 95 

12 CONTINUE B 96 

DO 13 1=1, N . B 97 

YOVSm=Y(n/Y(NI B 98 

13 CONTINUE B 99 

OELBA=Y(NI*RB B 100 

GO TO 15 B 101 

14 CONTINUE B 102 

CALL UPDATE B 103 

15 CONTINUE B 104 

WRITE (6,211 6 105 

CALL PRINTIT B. 106 

WRITE (6,221 B 107 

WRITE (6,23) CELTA,OELBA,REY,EMUSHOK B 108 

WRITE (6,241 B 109 

DO 16 1=1, N B 110 

WRITE (6,25) ETA( I ) ,ROV( I ) ,A( I ) ,P( I ),RO(I ),M( I ),T(I ) , AFOR( I ) B 111 

16 CONTINUE B 112 

RETURN B 113 

C B 114 

C B 115 

17 FORMAT (//8X,12HBOOY R AD IUS= , El 3 .5 , 25H CM, FREESTREAM VEL0CITY=,E1 B 116 

13.5, 17H CM/SEC, 0ENSITY=,E13. 5.5H G/CC/9X,18HA SHOCK (NON DIM)=,F8 B 117 

2.5, 19H, BETA (NON DIM ) =, F8. 5 , IH, , 16X, 9HPRESSURE- , E13 .5, 12H DYNES B 118 

3/SQ CM/9X,24HA8LATI0N RATE (NON 01 M ) = ,F8. 5 , BOH, MASS FRACTIONS- B 119 

4HYDRGGEN =,F7.5,10H, OXYGEN= ,F7 .5, 13H, NITROGEN =,F7.5,11H, CAR B 120 

5B0N =,F7.5) B 121 

18 FORMAT (9X,28HFREESTREAM COMPOSITION- N= ,F7.5,7H, 0= ,F7.5,7H, B 122 

1 C= ,F7.5,7F, H= ,F7.5) B 123 

19 FORMAT (50X,14HHITH RADIATION) B 124 

20 FORMAT I IHl , ///38X, 38HV I SCOUS STAGNATION STREAMLINE SOLUTION/) B 125 

21 FORMAT (1H1,////11X,40HINITIAL EQUILIBRIUM COMPOSITION PROFILES/) B 126 

22 FORMAT (1H1,////11X,21HINITIAL FLOW PROFILES) B 127 

23 FORMAT (//5X,6HCELTA=E16.8,2X,6H0ELeA=E16.8,2X,4HREY=E16.8,2X,21H B 128 

1 VISCOSITY(SHGCK)= ,E16.8) B 129 

24 FORMAT (4X,6HETA( I ) ,8X,6HR0V( I ) ,9X,4HA( I) ,10X,4HP( I ) ,10X,5HR0(I ),8 B130 

1X,4HH( I ) ,'8X,4HT( I )8X,7HAF0R( I ) ) B 131 

25 FORMAT (8E14i6) B 132 

END B 133- 
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APPENDIX A - Continued 


SU0ROUTIN6 RANHUG C I 

COHMGN /CAMFI/ AMFC (100 ) t AMFN( 100 I , AHFQ 1 100 ) t AHFH( 100 ) C 2 

COMMON /RANK/ RCSHOCK , VSHOCK , HSHOCK ,T SHOCK C 3 

COMMON /MISCL/ N,IOIM,CURV,DETAtOELTA,PSHOCK.REY,HW,EP,EROVfEA,MAX C t, 

ITIME.MLDIFPT.NICtERO C 5 

COMMON /FRSTRM/ H IN,P I N, RHOI NP.U INP , R8 C 6 

COMMON /MFINF/ ALP I NFC, ALP INFN, ALP I NFOt AtPINFH C 7 

COMMON /MFINJ/ AL P IN JC, ALP I N JN, ALP I NJC, ALP IN JH C 8 

COMMON /REFVALU/ EMUSHOK C 9 

EXTERNAL FOFXRO C 10 

EXTERNAL VISCCS C 11 

PSH0CK=0.95 C 12 

VSH0CK=1.+PIN-PSH0CK C 13 

HSHOCK=HIN-( VSH0CK*VSH0CK/2. I C lA 

CALL DIMPROP (PSHOCK, HSHOCK, PSHOKP.HSHOKPI C 15 

AINFC=ALPINFC C 16 

AINFN=ALPINFN C 17 

A1NF0=ALPINF0 C 18 

TSH0KP=14000. C 19 

TShOCK=TSHOKP C 20 

ROSHCCK=RHGSV(PSHOKP,HSHCKP*3.A83E-07,AINFC,AINFN,AINFO,TSHOKPI C 21 

CALL ITRl (R0SHCCK,.1,F0FXRG, .001, .001,20, ICODEI C 22 

IF (ICODE.EQ.OI GO TO 1 C 23 

IF (ICOOE.EQ.lt WRITE (6,21 C 24 

IF (IC0DE.E0.2i WRITE (6,31 C 25 

STOP C 26 

1 WRITE (6,41 PSHCCK,R0SH0CK, VSHOCK, HSHOCK C 27 

C C 28 

EMUSHOK=VISCOS(TSHOKP,PSHOKP» C 29 

RETURN C 30 

C 31 
C 32 

FORMAT (5X,37HMAXIMUM ITERATIONS EXCEEDED IN RANHUG) C 33 

FORMAT (5X,25H2ER0 DERIVATIVE IN RANHUGI C 34 

FORMAT (////5X,15HRANHUG SOLUTION//12H PSH0CK=E16 .8, 2X, 8HR0SH0 C 35 
lCK=E16.e,2X,7HVSH0CK=E16.8,2X,7HHSH0CK=E16.8) C 36 

ENC C 37- 


FUNCTION FOFXRO (ROSHOCI 0 1 

COMMON /RANK/ RQ SHOCK , VSHOCK , HSHOCK ,T SHOCK D 2 

COMMON /FRSTRM/ HIN, P IN, RHOI NP , UINP ,RB 0 3 

COMMON /MFINF/ ALP I NFC , ALP INFN, ALP I NFO, ALPINFH 0 4 

COMMON /MISCL/ N, IDIM,CURV,DETA, DELTA, PSHOCK, REY,HW,EP,EROV,EA, MAX D 5 

1TIME,ML0IFPT,NIC,ER0 0 6 

VSHOCK=1./ROSHOC D 7 

PSH0CK=PIN*1. -VSHOCK D 8 

HSH0CK=HIN-(VSHCCK*VSH0CK/2.) D 9 

CALL DIMPROP (PSHOCK, HSHGCK,PSHOKP,HSHOKP» 0 10 

AINFC=ALPINFC D 11 

AINFN=ALPINFN 0 12 

AINFO=ALPINFO 0 13 

TSHOKP=TSHOCK 0 14 

FOFXRO=RHOSY (PSH0KP,HSH0KP*3.483E-07, AINFC, AINFN,AINFC,TSHOKP) 0 15 

TSHOCK=TSHOKP D 16 

RETURN D 17 

END 0 18- 
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APPENDIX A - Continued 


SUBROUTINE UPDATE E 1 

COMMON /CA/ AdODtBdOll E 2 

COMMON /CY/ Y(101»,Y0YSd011 E 3 

COMMON /CETA/ ETAdOU E 4 

COMMON /CRO/ ROdOn E 5 

COMMON /CV/ VdOl) E 6 

COMMON /CP/ PdOlI E 7 

COMMON /CROV/ ROVdOK E 8 

COMMON /CH/ HdOU E 9 

COMMON /CSCH/ SCHdOll tSCl dOll E 10 

COMMON /CTEMP/ TEMP (101 1 , T IMP dOl » E 11 

COMMON /CT/ TdOll E 12 

COMMON /CAMFI/ AMFC (100 I , AMFN( 100 > , AMFOCIOO I , AMFHI 100 ) E 13 

COMMON /CFORONC/ AFOR ( 201 1 . AQNC( 2011 E 14 

COMMON /FRSTRM/ H IN,P I N, RHQI NP,U INP, R B E 15 

COMMON /MFINF/ ALP INFC t ALP INFN, A1.P INFO, ALPINFH E 16 

COMMON /MFINJ/ ALP I N JC , A LP I N JN, ALP IN JO, AL P 1 N JH E 17 

COMMON /MISCL/ N,IOIM,CURV, CETA, DELTA, PSHOCK,REY,HW,EP,EROV,EA, MAX E 18 

1TIME,MLDIFPT,NIC,ER0 E 19 

COMMON /RADCOMP/ DELBA,IRFLUX E 20 

NAMELIST /NAMl/ RO,ROV E 21 

NAMELIST /NAM2/ A,H,P,DELTA E 22 

namelist /NAM3/ AFOR E 23 

READ (5, NAMl) E 24 

READ (5,NAM2I E 25 

READ (5,NAM3) E 26 

DO 1 1=1, N E 27 

1 TEMPI 1 )=1./R0( I ) E 28 

CALL TINT (N,DETA,TEMP,TIMP) E 29 

DO 2 I=1,N E 30 

Y( I l=OELTA*TIMP( I) E 31 

SCH( n = l.*CURV'*OELTA*TIMP( I) E 32 

SCI ( I ) = SCH( I ) E 33 

2 CONTINUE E 34 

00 3 1=1, N E 35 

YOYS( I )=Y( n/Y(N) E 36 

3 CONTINUE E 37 

DELBA=Y(NI*RB E 38 

DO 4 1=1, N E 39 

4 AONC( I )=1.-AFCR( I) E 40 

DO 5 1=1, N E 41 

AMFC ( I) =AFOR ( I ) AALPINJC + AONCI I )*ALP1NFC E 42 

AMFN ( I )=AFOR( I l*ALPINJN*AONCd )*ALPINFN E 43 

AMFO(I)=AFOR( I)*ALPINJO+AONC(I)*ALPINFC E 44 

AMFH(I)=AFOR( II*ALPINJH+AONC(I)*ALPINFH E 45 

5 CONTINUE E 46 

RETURN E 47 

ENC E 48- 
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APPENDIX A - Continued 


SUBROUTINE FLOW F 1 

COMMON /CA/ AUOU.BdOl) F 2 

COMMON /CBETA/ BETAdOll F 3 

COMMON /CY/ YdOUtYOYSIlOU F 4 

COMMON /CRO/ ROdOl) F 5 

COMMON /CETA/ ETAdOll F 6 

COMMON /CV/ VdOn F 7 

COMMON /CP/ PdOl) F 8 

COMMON /CROV/ ROVdOU F 9 

COMMON /CDV/ CVdOl) F 10 

COMMON /CSCH/ SCMdOl I ,SCI (101 1 F 11 

COMMON /CTEMP/ TEMP (101 ) t TIMPdOl ) F 12 

COMMON /CEMU/ EMUdOlt F 13 

COMMON /CROMU/ ROMUdOl) F 14 

COMMON /CCAPH/ CAPHdOl) F 15 

COMMON /CH/ HdOl) F 16 

COMMON /CPRAN/ PRANdOll F 17 

COMMON /COR/ ORdOl) F IB 

COMMON /CXRO/ XNROdOl) F 19 

COMMON /CXH/ XNFdOl) F 20 

COMMON /CT/ TdOlt F 21 

COMMON /CAMFI/ AMFC ( 100 I . AMFNdOO 1 1 AMFO ( 1001 1 AMFHdOO » F 22 

COMMON /CFOROKC/ AFOR (201 1 . AONC( 201 I F 23 

COMMON /MISCL/ N, IDIM,CURV,0ETA, DELTA, PSHOCK,REY.HW,EP,EROV,EA, MAX F '24 

1TIME,MLDIFPT,NIC,ER0 F 25 

COMMON /MFINF/ ALP INFC, ALP I NFN, ALP INFO, ALP INFH F 26 

COMMON /MFINJ/ ALP IN JC , ALP I N JN , A LP I NJQ , ALP IN JH F 27 

COMMON /RADCGMP/ OELBA,IRFLUX F 28 

COMMON /DAMP/ I ROOAMP, ROCAMP , HOAMP F 29 

COMMON /C30/ CURVZ,BS,BETAZ(10n F 30 

COMMON /REFVALL/ EMUSHOK F 31 

DIMENSION TAMPdOU F 32 

UM0M=0 F 33 

NTIMES=0 F 34 

1 CONTINUE F 35 

WRITE (6,8) F 36 

WRITE (6,9) DELTA, OELBA.REY F 37 

WRITE (6,10) F 38 

DO 2 1=1, N F 39 

WRITE (6,11) ETA(I),V(I),QR(I),CAPH(I),V(I),EMU(I),PRAN(I),YCYS(1) F 40 

2 CONTINUE F 41 

WRITE (6,12) F 42 

DC 3 .1 = 1, N F 43 

WRITE (6,13) ETA(I),ROV(I),A(I),P(I),RO(I),XNRCd),H(I),XNH(l) F 44 

3 CONTINUE F 45 

WRITE (6,14) F 46 

DO 4 1=1, N F 47 

WRITE (6,13) ETA( I ),T( I ) ,AMFC( I ) ,AMFN( I ),AMF0(1I,AMFH( I ) ,AFOR( I ) ,8 F 48 

1(1) F 49 

4 CONTINUE F 50 

CALL SET (N,RCV,TAMP) F 51 

CALL MASS F 52 

ICCNT=ITEST(N, BOV, TAMP, EROV) F 53 

CALL SET (N,P,TAMP) F 54 

CALL YMOMNTM F 55 

IYM0M=ITEST(N,P,TAMP,EP) F 56 

CALL SPECIES F 57 

CALL VISPRAN F 58 

CALL SET (N, A, TAMP) F 59 

CALL XMCMNTM ( A , BET A, CUR V ) F 60 

IXMOM=ITEST(N,A,TAMP,EA) F 61 

IF ( IDIM.NE.3) GO TO 5 F 62 

CALL SET (N, 8, TAMP) F 63 

CALL XMCMNTM I A , BET AZ, CURVZ ) F 64 

IZMOM=ITEST(N,B,TAMP,EA) F 65 

5 CONTINUE F 66 

CALL ENERGY F 67 

IF ( dCCNT.NE.O).OR.(IYMGM.NE.O) .OR.dXMCM.NE.O).OR.( IZMOM.NE.O) ) F 68 

IGO TO 6 F 69 

GO TO 7 F 70 

NTIMES=NTIMES+1 F 71 

IF (NTIMES.LE.MAXTIME) GC TO 1 F 72 
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WRITE (6fl6) f 73 

STOP F 74 

7 CONTINUE F 75 

WRITE 16,151 F 76 

RETURN F 77 

C F 78 

C F 79 

8 FORMAT (1H1,////11X,13HFL0W PROFILES) F 80 

9 FORMAT (///5X,6H0ELTA=E16.8,2X,6H0ELBA=E16.8,2X,4HREY=E16.8/I F 81 

10 FORMAT (4X,6HETA( n,9X,4HY(II ,1CX,5HCR( 1 1, 8X,7HCAPH( I ) ,10X,4HV( I) , F 82 

19X,6HEMU( I 1 ,8X,7HPRAN(I ) ,8X,7HYOYS( I) ) F 83 

11 FORMAT (8E14.6) F 84 

12 FORMAT (4X,6HETA( I ) ,8X,6HR0VI I ) ,9X,4HA( I) ,10X,4HP( I ) ,10X,5HR0( I ) ,8 F 85 

lX,7HXNROm ,9X,4HH( I),9X,6HXNH(I n F 86 

13 FORMAT (8E14.6) F 87 

14 FORMAT (4X,6HETA( I) ,9X,4HT( I) ,9X,7HAMFC(I I ,8X,7HAMFN( I ),8X,7HAMF0( F 88 

in ,6X,7HAMFH( I ) ,7X,7HAFOR( I) ,7X,4HB( I ) ) F 89 

15 FORMAT (5X,23HFLOW SOLUTION CONVERGED) F 90 

16 FORMAT (5X,27HFLOW SOLUTION NOT CONVERGED) F 91 

END F 92- 


SUBROUTINE CHEXRO nCHXRO) G 1 

COMMON /CH/ HdOl ) G 2 

COMMON /CP/ P(lOl) G 3 

COMMON /CRO/ RC(lOl) G 4 

COMMON /CXRO/ XNRO(lOl) G 5 

COMMON /CAMFI/ AMFC ( 100 ) , AMFNI 100 ) , AMFOI 100 ) , AMFHI 100 ) G 6 

COMMON /MISCL/ N , I DIM ,CURV , OET A , CELT A, PSHOCK , REY, HW, EP , EROV, EA, MAX 6 7 

1TIME,MLDIFPT,NIC,ER0 G 8 

COMMON /DAMP/ IROOAMP.RODAMP.HOAMP G 9 

DIMENSION ABRAT(lOl), RATIO(lOl) G 10 

DIMENSION PP(lOl), HP(lOl) G 11 

DIMENSION AMFCPIIOO), AMFNP(IOO), AMFOPIlOO) G 12 

EXTERNAL RHOSY G 13 

ICHXR0=0 G 14 

CALL SET (N,AMFC,AMFCP) G 15 

CALL SET (N,AMFN, AMFNP) G 16 

CALL SET (N,AMFO,AMFOP) G 17 

TGUES=4000. G 18 

DO 1 1=1, N G 19 

CALL DIMPROP (Pm,HlI),PPm,HPm) G 20 

XNRO(I)=RHOSY(PP(I),HPm*3.483E-07,AMFCP(I),AMFNP(I),AMFOP(I),TGU G 21 

lES) G 22 

1 CONTINUE G 23 

DO 2 1=1, N G 24 

RATIO(I)=(XNRCm-RO(I))/XNROm G 25 

ABRATd )=ABS(RATIO( I ) ) G 26 

IF ( ABRATd ) .GT.ERO) ICHXR0=1 G 27 

2 CONTINUE G 28 

IF dPODAMP.GT.O ) GO TO 4 G 29 

CO 3 1=1, N G 30 

3 ROd )=XNRO(I ) G 31 

IR0DAMP=1 G 32 

GO TO 6 G 33 

4 CONTINUE G 34 

DO 5 1=1, N G 35 

ROd )=RCOAMP*RO( I )♦( l.-RCDAMP)*XNRO(I ) G 36 

5 CONTINUE G 37 

6 CONTINUE G 38 

RETURN G 39 

END G 40- 
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c 

c 

c 

c 


1 

2 

3 

4 

5 

6 


7 

8 

9 

10 
11 


12 

C 

C 

C 

13 

14 


15 


16 


17 

18 

19 

20 
21 
22 


SUBROUTINE MASS 

COMMON /CA/ A(lOl) ,8(101 » 

COMMON /CETA/ ETA(lOl) 
common /CY/ Y(101» ,Y0YS( 1011 
COMMON /CRO/ RO(lOl) 

COMMON /CROV/ ROVClOlt 

COMMON /CSCH/ SCH ( 101 ) , SC 1 1 101 ) 

COMMON /CTEMP/ TEMP ( 10 1 1 , T IMP (101 1 
COMMON /CV/ VdOll 

COMMON /MISCL/ N, IOIM,CURV,OETA,OELTA,PSHOCK,REY,HW,EP,EROV,EA,MAX 
1TIME,ML0IFPT,MC,ER0 
COMMON /FRSTRM/ H I N , P I N , RHOI NP , U I NP , R B 
COMMON /RADCOMP/ DELBA.IRFLUX 
COMMON /C30/ CURVZ.BS.BETAZdOl I 

MASS CONSERVATION PACKAGE 
COMPRESSIBLE ETA TRANSFORMATION 

IF ( IDIM.EQ.2I GO TO 2 
IF ( I0IM.EQ.3 1 GO TO 4 
DO 1 1=1, N 
TEMPI! I=A( I) 

GO TO 6 
DO 3 1=1, N 

TEMPI I )=2.*A( n*SCH( I 1 
GO TO 6 
CO 5 1 = 1, N 

TEMPIl 1=A( II*SCH( 1 )+B( n*SCl ( I 1 
CONTINUE 

CALL TINT IN, CETA, TEMP, TIHP) 

IF ( IDIM.E0.2) GO TO 7 
IF ( IDIM.EQ.31 GO TO 9 
CALL SET (N, SCH, TEMPI 
GO TO 11 
DO 8 1=1, N 
TEMPI I )=SCH( I )**2 
GO TO 11 
DO 10 1=1, N 
TEMPII)=SCHI IlMSCIin 
CONTINUE 

DELTA=I TEMPI NI-ROVIl I l/TIMPINI 
DO 12 1=1, N 

BOV I I t = IOELTA*TIMPI I 1 +ROV( 1 1 I /TEMP I 1 1 
DO 13 1=1, N 

MASS CONSERVATION SATISFIED 

VI n=Rovi 1 1/Rci n 

DO 14 1=1, N 

TEMPI I 1=1. /ROI II 

CALL TINT (N,DETA,TEMP,TIMPI 

00 15 1=1, N 

Y(I)=0ELTA*TIMP(I I 

SCH I I l = l.+CURV*OELTA*TIMP( II 

CONTINUE 

DO 16 1=1, N 

Y0YSII)=YII)/YINI 

CONTINUE 

DELBA=Y(N)*R8 

IF I ICIM.EQ.2) GO TO 18 

IF I IDIM.EQ.31 GO TO 20 

00 17 1 = 1, N 

SCI I I ) = 1. 

GO TO 22 
DO 19 1=1, N 

scim=scHin 

GO TO 22 
DO 21 1 = 1, N 

SCI! I ) = 1.+CURVZM0ELTA*TIMPI I I 

CONTINUE 

RETURN 

END 
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SUBROUTINE YMGMMM II 

CORMCN /CETA/ ETA(lOl) I 2 

COMMON /CP/ PdOll I 3 

COMMON /CRO/ ROdOlJ I 4 

COMMON /CRCV/ RCVdOll I 5 

COMMON /CDV/ CVUOH I 6 

COMMON /CTEMP/ TEMP ( 10 1 ) f OP dO 1 ) I 7 

COMMON /CV/ VdOl) I 8 

COMMON /MISCL/ N , I DI M , CURV , OET A.OELT A .PSHQCK.REY, HW, EP. EROVt EA ,MAX I 9 
lTIME,MLOIFPT,NICtERC I 10 

THIS SUBROUTINE INTEGRATES THE Y MOMENTUM EQUATION TO GIVE I 11 

P(YI I 12 

CALL DIET (V.ETA.DVtN) I 13 

DO 1 I=ltN , I 14 

1 TEMPd )=ROV( I )*DV( I) I 15 

CALL TINT (N.DETA, TEMP, DPI I 16 

PSTAG=PSHOCK*DP«N) I 17 

DO 2 1=1, N I 18 

P( I l = PSTAG-OP( II I 19 

2 CONTINUE I 20 

RETURN I 21 

END I 22- 
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SUBRCUTINE SPECIES J 

COMMON /CRO/ ROaOl) J 

COMMON /CT/ TdOll J 

COMMON /CETA/ ETA(lOl) J 

COMMON /CROV/ ROVIIOU J 

COMMON /CSCH/ SCHIlODtSCKlOl) J 

COMMON /CD12/ 012(101) J 

COMMON /CAMFI/ AMFC ( 100 ) , AMEN ( 100 I t AM FO ( 100 1 1 AMFH( 100 ) J 

COMMON /CFORONC/ AFOR ( 201 ) , AONC ( 201 1 J 

COMMON /MISCL/ N , 1 01 H , CURV , DET A , CELT A , PSHOCK , REY , HW t EP , EROV , EA , MAX J 
ITlMEtMLDIFPT.NIC.ERO J 

COMMON /MFINJ/ AL P 1 N JC . ALP I N JN, ALP I NJQ , AL P 1 NJH J 

COMMON /MFINF/ AL P 1 NFC , ALP INFN , ALP I NFO, ALP 1 NFH J 

DIMENSION AS(lOl). BS(lOl), CS(lOl), CS(lOl) J 

DIMENSION RK(201). 012K(201). RVK(2C1I. SHK(201) J 

DIMENSION ETADIF(201I. AK(201). BK(201lt CK(201). DK(201) J 

DIMENSION DK201), 00(201) J 

CALL CIFCOEF J 

AINF=1.0 J 

AINJ=1.0 J 

NLI = N-:1 J 

CO 1 1 = 1 , N J 

DI ( I )=SCH( n*RO( I 1*R0(I )*D12( I l*SCI( 1 I J 

C0( I )=DETA*SCH( I )*ROV( I )*SCI ( 1 ) J 

BOUNDARY CONDITIONS - J 

DIFFUSION OF SHOCK LAYER SPECIES TO WALL PERMITTED J 

AS(1)=0. ' J 

BS(1 )=-DI(l )»C0(1 )*DELTA J 

CS(1)=0I(1I J 

DS(1 )=DELTA*CO( 1 )*AINJ J 

AS(NI=0. J 

BS(N)=1.0 J 

CS(N)=0. J 

0S(N)=0 J 

J 

00 2 1=2, NLI J 

SIGN=0. J 

IF (ROVm.GT.O.I SIGN=1.0 J 

AS(I)=DI( I-1)*DI( I )-2.*CC( I)*DELTA*(1.-SIGN) J 

BS(I)=-(01 ( I-1)»2’*0I( I)+OI ( dl) ) J 

BS(I ) = 8S( I l*2.*C0( I )*DELTA*( 1. -SIGN 1-2. *C0( II*OELTA*SIGN J 

CS( I )=0I( I (♦■OK I+l l♦2.♦C0(n*DELTA*SIGN s J 

DS(I)=0. J 

CALL TRIOIAG ( N, AS , BS , CS ,0S I J 

CO 3 1=1, N J 

AFCRm=DS(I) J 

AONC ( n=l .-AFOR( I 1 J 

CO 7 1=1, N J 

IF ( AFOR( I I .LT.l.E-06 I A, 5 J 

WRITE (6,10) J 

WRITE (6,11) ETA( I ) , AFOR( I ) J 

AFCR(I)=0. J 

A0NC(I)=1.0 J 

CONTINUE J 

IF ( AFOR( I ) .GT..S99999) 6,7 J 

WRITE (6,121 J 

WRITE (6,11) ETA( I ) ,AFOR( I ) J 

AFCR(I)=1.0 J 

A0NC(I)=0. J 

WRITE (6,11) J 

CONTINUE J 

DO 8 1=1, N J 

AMFC( I )=AFOR( I )*ALPINJC+AONC( I )*ALPINFC J 

AMFNI I)=AF0R(I)*ALP1NJN«-A0NC(I l^ALPINFN J 

AMFO( I )=AFOR( I ) ♦ ALP I N J0+ AONC ( I )*ALPINFO J 

AMFH( I )=AFOR( I ) *ALP I NJH4-A0NC ( I )♦ ALP I NFH J 

AMFC( I )=AMFC ( I )/(AMFC( I) ♦AMFN( I )*AMFO( I l+AMFHl ID J 

AMFN( I )=AMFN( I )/(AMFC( I ) tAMENI I )*AMF0( I )*AMFH( II) J 

AMFO( I )=AMFO( I )/ ( AMFC( I )+AMFN( I )+AMFO( I ) + AMFH( I ) ) J 

AMFH( 1 ) = AMFH( I 1/ ( AMFC( I ) *AMFN( I I tAMFO( I )*AHFH( I ) ) J 

WRITE (6,131 J 
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DO 9 1=1, N J 72 

9 WRITE (6,1A) ETAm,AMFCm,AMFN(n,At<FO<n,AMFHm J 73 

RETURN . _ J 74 

J 75 
J 76 

10 FORMAT (1H0,5X,24HAF0R( I) LESS THAN l.E-06) J 77 

11 FORMAT (2E14.6J J 78 

12 FORMAT ( IHO, 5X,28HAFOR( I) GREATER THAN .999999) J 79 

13 FORMAT (//10X,18HELEMENTAL PROF I LES/5X ,6HETA(I ) , 6X , 9HC ARBONI I) , 4X , J 80 

lllHNITRCGENI 1 1 ,3X,9HOXYGEN( n ,5X,11HHY0R0GEM I)/) J 81 

14 FORMAT (8E14.6) J 82 

END J 83- 
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SUBROUTINE XMORNTK (A.BETA.QI 
COMRON /CETA/ ETAdOll 
COMMON /CRO/ ROtlOll 
COMMON /CROV/ RCVdOll 
COMMON /CDV/ CVdOU 
COMMON /CEMU/ EMUdOl) 

COMMON /CROMU/ ROMUdOl) 

COMMON /CSCH/ SCHdOllfSCIdOlI 

COMMON /MISCL/ N,IDIM,CURV,OETA,OELTA,PSHOCK,REY.HW,EP,EROV,EA,MAX 
ITIME.MLDlFPT.NICtERO 

THIS SUBROUTINE SOLVES THE VISCID X MOMENTUM EQUATION 

DIMENSION AMdODt BMdODt CHdOlIt CMdDDt AdOlIt BETAdOlt* T 
lUMPdOl), TIMPdOlIt AVMdOl), BVM(lOl), CVMdOlK DVMdOll 
DIMENSION DELdOll 
TEMP=l./(2.*OELTA*DETA» 

NL1=N-1 

DELT=l./(2.*REY*OELTA*OETA*OETAI 

Ad(=0. 

NTIMES=0 
NXIT= 10 
EPSX=.01 
CONTINUE 

IF (NTIMES.LT.NXIT I GO TC 2 
WRITE (6,8) 

STOP 

CONTINUE 

NTIMES=NTIMES*1. 

AMd)=0. 

AM(NI=0. 

6Md)=l. 

BM(N)=1. 

CMd)=0. 

CM(N)=0. 

DMd)=0. 

0M(NI=0. 

AVMd)=0. 

AVM(N)=0. 

8VM( 1»=1. 

BVM(NI=1. 

CVMdl=0. 

CVM(N)=0. 

DVMd)=0. 

DVM(NI=0. 

TUMPd)=0. 

TIMPd)=0.0 
DO 3 1=2, NLl 

CM (I l=-ROV( I l♦TEMP♦SCId I^ROd ) 

AMd l=-CM( I ) 

CONTINUE 
DO 4 1=1, N 
TUMP( I) = 1./R0(I ) 

CALL TINT (N,DETA,TUMP,TIMP) 

DO 5 1=1, N 

ROMU( I )=R0( I l*EMU( I I 
DEL( n = DELT*SCI( I l*RO( I ) 

DO 6 1=2, NLl 

BM( I ) = 2.*R0( I )*A( I )-RCV( I 1*0 

CM< I )=8ETA( I )-R0( I )*A( I )**2*R0V( I l*(A d*l)-A( 1-1 ))*TEMP*SCI ( I (♦R0( 
in4R0V(I)*Am*Q 

AVM(I ) = -(ROMU( I-1)*R0MU( I) )*DEL( I) /DELTA 

BVM( I )=(ROMU( I-l |42.*R0MU( I ) +ROMU ( I ♦ 1 ) ) ♦DEL ( I) /DELTA 

CVMd )=-(ROMU( 1 )+ROMU( I + ll )*DELd ) /DELTA 

DVM( I ) = -AVM( I)#A( I-1I-BVM( I ) «A (I I -C VM d ) *A ( I *1 ) 

AMd ) = AM( I )*AVM( 1 ) 

BMm=BMm+BVMm 
CM( I )=CM{ I )4CVM( I) 

0Mm=DM(I)4DVMm 

CONTINUE 
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CALL TRIDIAG ( N , AM, BM, CM ,DH ) K 71 

IG0R=0 K 72 

DO 7 1=1, N K 73 

A( I l=A( I )*0M( n K 7A 

IF ( ABSIDMI I ) I.LT.EPSX) GO TO 7 K 75 

IG0R=IGCR*1 K 76 

7 CONTINUE K 77 

IF (IGOR.NE.O) GO TO 1 K 78 

RETURN K 79 

K 80 
K 81 

FORMAT (10X,28HX MOMENTUM FAILS TO CONVERGE! K 82 

ENC K 83- 
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SUBROUTINE ENERGY 

WINDWARD DIFFERENCE FORM OF ENERGY ECUATION 

ROmiS UPDATED IN THIS SUBROUTINE BEFORE H IS UPDATED 
COMMON /CA/ AIIOII.BIIOI) 

COMMON /CY/ Y(101),YOYS(101I 
COMMON /CRO/ RO(lOl) 

COMMON /CETA/ ETAdOll 
COMMON /CV/ VIlOll 
COMMON /CP/ P(lOl) 

COMMON /CT/ TdOlI 
COMMON /CROV/ ROV(lOl) 

COMMON /COV/ DV(lOl) 

COMMON /CSCH/ SCHIIOD.SCKIOI) 
common /CTEMP/ TEMP(101»,TIMP(101I 
COMMON /CEMU/ EMUdOll 
COMMON /CROMU/ ROMUdOl) 

COMMON /CPRAN/ PRANdOl) 

COMMON /CCAPH/ CAPH(lOl) 

COMMON /CH/ HdOlI 
COMMON /CQR/ CRdOX) 

COMMON /CDOR/ OQRdOll 

COMMON /CXH/ XNHdOll 

COMMON /RADCOMP/ OELBA.IRFLUX 

COMMON /FRSTRM/ H I N , P I N, RHQ I NP, U INP , R fi 

COMMON /MISCL/ Nf lOIM.CURV.OETAf OELTA.PSHOCK ,REY,HW,EP,EROV,EA,MAX 
ITIME.MLDIFPT.NIC.ERQ 
COMMON /DAMP/ IROOAMP.ROOAMP.HOAMP 
DIMENSION AMdOlIf BMdOlIt CMdOll, CMdOll 
DIMENSION VISdOlIf CQNVdOlI, VGdODt SCHSdOU 
DIMENSION OVMdOl), OQRMdOll 
CALL RAOARAY 
DO 1 1 = 1, N 

RQMU( I l=RO( I >*EMU( I I 

CONTINUE 

00 2 1=1, N 

scHS(ii=scHm*sci(n 

VIS( I )=SCHS( I »*ROMU(I l/(PRAN(II*OETA) 

CONV( I l=DELTA*REY*SCHS( I )*ROV( I » 
vG( n=vis( I i*Rov( n/ROd > 

CONTINUE 
NLI=N-1 
DO 3 1=2, NLI 
SIGN=0. 

IF (ROV( n .GT.O. I SIGN = 1.0 

AMd »=VIS( I-H+VISd l-2.*C0NV( I )*d.-SIGN) 

BMd )=-( VISI I-H »2.*VI S( I) ♦VISd»m 

BM( I ) = BMd l+2.*C0NV( I )*d.-SIGNI-2.*CQNV( I »*SIGN 

CMd l=VISd •♦VISd*l) + 2.*C0NV( DESIGN 

DVM( I) = (VG( I-l dVGtl ) )*Vd-ll-(VGd-n*2.*WGm*VGt I ♦ 1 d d I ♦( VG( 
II ) »V6( dll t «V( dll 

OQRMd )=OELTA*REY*(-SCHS( I I *QR d -1 d ( -SCHSl I-1)»SCHS( dl d*ORd dS 
ICHSI IdORd+ll ) 

DMd 1=DVM( I l+DORMI I ) 

CONTINUE 
AMd 1=0. 

AM(NI=0. 

BMd)=l. 

BM(NI=1. 

CMddO. 

CM(N)=0. 

DM(1 |=HW<-Vd dVd )/2. 

CM(NI=HIN 

CALL TRIDIAG ( N , AM , BM, CM, DM I 

CO 4 1 = 1, N 

CAPHddDMm 

XNH( I l = CAPH( I )-Vd l*Vd 1/2. 

H( I l=HOAMP*H( I dCl.-HOAMPdXNHd I 
CONTINUE 
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DO 6 I=lfN L 71 

IF (Hm.LT.HUn 5,6 L 72 

5 CONTINUE L 73 

WRITE (6,71 L 7A 

WRITE (6,8» ET/m,Hm 1. 75 

H(I)=Hm L 76 

6 CONTINUE L 77 

RETURN ' L 78 

L 79 
L 80 

, FORMAT (1H0,5X,19HH( I) LESS THAN H(ll» L 81 

FORMAT (2E1A.6) L 82 

END L 83- 
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SUBROUTINE RADARAY 

RADARAY GIVES THE RADIATION HEAT ELUX AT NIC POINTS 
ANC LINEARLY INTERPOLATES FOR THE REMAINING N POINTS 


N MUST BE LESS THAN 100 

COMMON /CRO/ ROClOll 
COMMON /CETA/ ETA(lOl) 

COMMON /CY/ YdOl liYOYSdOU 
COMMON /CP/ PdOll 
COMMON /CH/ HdOU 
COMMON /COR/ QRdOlI 
COMMON /CXRO/ XNROdOlI 
COMMON /CT/ TdOll 

COMMON /CAMFI/ AMFC dOO I , AMFN dOO » t AMFOdOO 1 1 AMFHdOO > 

COMMON /RADCOMP/ OELBA.IRFLUX 

COMMON /MISCL/ N, I DI M, CURV , OETA, OELT A , PSHOCK.RE Y, HW, EP, EROV, EA, MAX 
ITIMEfMLCIFPT.MCtERO 
COMMON /FRSTRM/ H I N , P I N, RHO I NP , U I NP , R 8 
COMMON /DAMP/ I ROOAMP ,RQOAMP , HOAMP 

DIMENSION PHdOOIf HHdOOl, ETAHdOOIt PHPdOOd HHPdOOIt ORHPdO 
10), RHOHP(IOO), NICNdOOl 
DIMENSION ETANdOOl, QRPdOOl 
NN=N 

DO 1 1=1. NN 
ETAHd )=Y( n/Y(NI 
PHd)=Pd) 

HHd ) = H( I ) 

CALL DIMPROP (PH ( I I ,HHd ) , PHP < I » , HHP ( I ) I 
HHPt I ) = HHP( !)♦( .2389E-07 ) 

NRA0=(NN-1)/(NIC-1I 
DO 2 1=1, NIC 
J = NRAD*(I-1»H 
ETAN( n = ETAH( J) 

CONTINUE 

DO 3 1=1, NIC 

NICN( I l=NRAO*( I-l )♦! 

IF ( IRFLUX.LE.O) GO TO 6 

CALL RATRAP ( CE L BA , N I C , N IC N, NN, 4000 .0 , E T AH.PHP , AMF C , AMFN, AMFO, HHP , 
1QRHP,1,2,RH0HP,1,T) 

DO 4 1=1, NN 

CALL FTLUP (ETAHd), ORPdl,l, NIC, ETAN.QRHPI 

CONTINUE 

DO 5 1=1, N 

OR (I ) = ORP( I )/(RH0INP*UINP*UINP*UINP)*l.E*07 

GO TO 8 

CONTINUE 

CALL RATRAP ( GEL8 A, N IC , N ICN, NN, 4000 .0 , ETAH, PHP , AMFC , AMFN, AMFC, HHP, 
1QRHP,0,2,RH0HP,1,T I 
CO 7 1=1, N 
QRd )=0. 

CONTINUE 

IF (IRODAMP.GT.O ) GO TO 10 
DO 9 1 = 1, N 

ROd )=RHOHP( I )/RHOINP 

GO TO 13 

CONTINUE 

DO 11 1=1, N - 

XNRO( I )=RHOHP( II/RHOINP 

DO 12 1=1, N 

ROd ) = R0DAMP4R0(I)-M1.-RCDAMP)*XNR0(I ) 

CONTINUE 

RETURN 

END 


M I 
M 2 
M 3 
M 4 
M 5 
M 6 
M 7 
M 8 
M 9 
M 10 
M 11 
M 12 
M 13 
H 14 
M 15 
H 16 
M 17 
M 18 
M 19 
M 20 
M 21 
M 22 
M 23 
M 24 
M 25 
M 26 
M 27 
M 28 
M 29 
M 30 
M 31 
M 32 
M 33 
M 34 
M 35 
M 36 
M 37 
M 38 
M 39 
M 40 
M 41 
M 42 
M 43 
M 44 
M 45 
M 46 
M 47 
M 48 
M 49 
M 50 
M 51 
M 52 
M 53 
M 54 
M 55 
M 56 
M 57 
M 58 
M 59 
M 60 
M 61 
M 62 
M 63 
M 64 
M 65 
M 66- 
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10 


SUBROUTINE PRINTIT 


N MUST BE LESS THAN 100 

COMMON /CRO/ RO(lOl) 

COMMON /CETA/ ETAdOlt 
COMMON /CY/ Y(101».Y0YS(101» 

COMMON /CP/ P(101» 

COMMON /CH/ HdOll 
COMMON /COR/ QRdOK 
COMMON /CXRO/ XNROdOl) 

COMMON /CT/ Td01» 

COMMON /CAMFI/ AMFC (100 ) » AMFNdOO ) , AMF0{ 100 I , AMFHdOO 1 
COMMON /RAOCOMP/ OELBA.IRFLUX 

COMMON /MISCL/ N, ID 1 M,CURV . OETA, DELTA ,PSHOCK,REY, HW , EP ,EROV, EA# MAX 
ITIMEf MLOIFPT.NICtERO 
COMMON /FRSTRM/ H IN , P IN, RHOI NP , UINP , RB 
COMMON /DAMP/ IRODAMP,ROCAMP,HOAMP 

DIMENSION PHdOO), HHdOOI, ETAHdOOI, PHPdOOl, HHPdOOI, QRHPdO 
101, RHOHPdOOI, NlCNdOOI 
DIMENSION ETANdOOl, QRPdOO) 

NN=N 

DO 1 1=1, NN 
ETAHd t=Y( n/Y(Nl 
PH(I )=PU ) 

HH( I l = H( I ( 

CALL DIMPROP (PH ( I) , HH ( I I , PHP( I ) , HHP ( I ) > 

HHP( I l=HHP( !)♦( .2389E-07 ) 

NRAD=(NN-1)/(NIC-1) 

DO 2 1=1, NIC 
J = NRAD*d-l)*l 
ETAN(H=ETAH( J) 

CONTINUE 
DO 3 1=1, NIC 
NICN(I)=NRA0*(I-1»+1 
IF (IRFLUX.LE.O) GO TO 6 

CALL RATRAP (DEL6A, NI C ,N ICN,NN, 4000 .0 , ETAH, PHP , AMFC, AMFN, AMFO, HHP , 
1QRHP,1,2,RHOHP,0,T> 

WRITE (6,9) 

DO 4 I=1,NN 

CALL FTLUP ( ET AH ( I ) , QRP ( II ,2 , NIC ,ET AN ,QRHP ) 

WRITE (6,101 ETA( 1 1 ,Y( I ) ,YOYS( I ) ,QRP( I I 

CONTINUE 

DO 5 1=1, N 

CR( I )=QRP( I )/ (RF0INP*UINP*UINP*UINP)*l.E+07 

GO TO 8 

CONTINUE 

CALL RATRAP ( DELBA , N I C ,N ICN, NN , 4000 .0 , ETAH , PHP , AMFC, AMFN, AMFO, HHP, 
1CRFP,0,2,RHOFP,0,T) 

CO 7 1=1, N 
QR(I)=0. 

CONTINUE 

RETURN 


FORMAT (4X,6HETA( I I ,9X,4HY( 1 1 ,8X,7HYCYS(I I , 3X , 15HQR AD ( I ),W/SQ CM) 

FORMAT (4E14.6I 

END 


N 1 
N 2 
N 3 
N 4 
N 5 
N 6 
N 7 
N 8 
N 9 
N 10 
N 11 
N 12 
N 13 
N 14 
N 15 
N 16 
N 17 
N 18 
N 19 
N 20 
N 21 
N 22 
N 23 
N 24 
N 25 
N 26 
N 27 
N 28 
N 29 
N 30 
N 31 
N 32 
N 33 
N 34 
N 35 
N 36 
N 37 
N 38 
N 39 
N 40 
N 41 
N 42 
N 43 
N 44 
N 45 
N 46 
N 47 
N 48 
N 49 
N 50 
N 51 
N 52 
N 53 
N 54 
N 55 
N 56 
N 57 
N 58 
N 59- 


FUNCTION D2DP2 (X) 0 1 

COMMON /CD20P2/ D2DP2C 0 2 

020P2=D20P2C 0 3 

RETURN 0 4 • 

END 0 5- 
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SUBROUTINE DIET (Y,X,OY,N) P 1 

DIMENSION YllOlIf X(lOl), DY(lOl) P 2 

DY(1) = (Y(2>-Ym)/(X(2I-X(l)) P 3 

DYtNI=lY(N)-YlN'in/IXlN)-X(N-l) J P ‘t 

NX=N-1 P 5 

CO 1 I=2fNl P 6 

X1=X(I-1)-X( I» P 7 

X3=X( I + ll-X( I ) P 8 

DY( I ) = ( (Y( I+l )-Y( I ) I*X1**2*( Yl n-Y( I-in*X3**2>/( (X1-X3»*X1*X3) P 9 

1 CONTINUE P 10 

RETURN P II 

END P 12- 


SUBROUTINE TRIOIAG (N,A,B,C,DI 
DIMENSION A(201l< B(20Ut C(201)t 0(231> 

THIS ROUTINE SOLVES AX=0 FOR A TRICIAGCNAL 

N1=N-1 

c(i)=cm/B(i) 

DJi)=Dm/Bm 
DO 1 M=2,N'l 

D(MI=(0(HI-AtMl*0(M-l( t/(B(M|-A(M|»C(M-l)| 

1 C(M|=C(M)/(8(M1-A(M1*C(M-1) ) 
0(NI=(D(N»-A(N»*D(N-in/(B(NI-A(NI*C(N-in 
DO 2 J=2.N 

M=N+1-J 

2 0(M|=D(MI-C(MI*0(M*1) 

RETURN 

END 


Q 1 
0 2 
0 3 

0 4 

Q 5 
0 6 
Q 7 
0 8 
Q 9 
0 10 
0 11 
0 12 
Q 13 
0 14 

0 15 

0 16 
Q 17- 


SUBROUTINE SET (N,A,B» R 1 

DIMENSION AdOUf B(lOl) R 2 

DO 1 1=1, N R 3 

1 Bm=A(i) R ^ 

RETURN R 5 

END 


FUNCTION ITEST (N,A,B,E) 

DIMENSION A(lOl), B(101» 

ITEST=0 
DO 1 J=1,N 

1 IF (ABSIAI J)-B( J) ) .GT.EI I TEST= I TEST* 1 

RETURN 
END 
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2 
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SUBRCUTINE DIFCOEF 
COMMON /CETA/ ETA(lOl) 

COMMON /CP/ P(lOl) 

COMMON /CH/ HdOll 
COMMON /CT/ TtlOl) 

COMMON /C012/ 012(1011 

COMMON /MISCL/ N, IOIM,CUPV,OETA,DELTA,PSHOCK,REY,HW,EP,EROV,EA,MAX 
ITl ME tMLOl FPT ,NIC f ERO 
COMMON /FRSTRM/ H I N , P I N , RHC 1 NP , U INP t R B 
COMMON /COINGll/ TSTR(18l,0MEGAll(l'8) 

CIMENSION TST(lOl), OMEGdJll 
CIMENSION PHPdOlIt HHPdOU 
REAL M1,M2,M12 

BINARY DIFFUSION COEFFICIENT 

1- ATOMIC NITROGEN 

2- ATOMIC HYDROGEN 
MI=1A.01 
M2=1.G08 
SIGMA1=3.298 
SIGMA2=2.708 
EPS0K1=71.4 
EPS0K2=37. J 

SIGMA12=. 5#(SIGMA1*SIGMA2» 

EPSOK12=SORT ( EPS0K1*EPS0K2 I 

DO 1 1 = 1, N 

TST( n = T( I )/EPS0K12 

CALL FTLUP ( T S T ( I) , OM.EG ( I I , -1 , 18 , TSTR , OMEGA H ) 

CONTINUE 
00 2 1=1, N 

CALL OIMPROP (P (I I ,H( I I , PHP( I I ,HHPd ) ) 

CONTINUE 

M12=(M1*M21/(2.mm1*M2I 

M12=S0RT(M12) 

CC=0.002628*M12/(SIGMA12*SIGMA12I 
DC 3 1=1, N 

012( n = OC/(PMP( I l*OMEGd I I *(T( I » J**1.5/(R8*UINP) 

CONTINUE 

RETURN 

ENO 


T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 

T 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 
39- 


SUEROUTINE VISPRAN U 1 

COMMON /CP/ PdOll U 2 

COMMON /CH/ HdOl) U 3 

COMMON /CT/ TdOll U 4 

COMMON /CEMU/ EMUdOll U 5 

COMMON /CPRAN/ PRANdOU U 6 

COMMON /CETA/ ETAdOll U 7 

COMMON /FRSTRM/ HIN ,PI N,RHOINP ,UINP ,R0 U 8 

COMMON /MISCL/ N, IOIM,CURV,OETA, DELTA, PSHOCK,REY,HW,EP,EROV,EA, MAX U 9 

lTIMe,MLDIFPT,NIC,ERO U 10 

COMMON /REFVALU/ EMUSHOK U 11 

DIMENSION PPdOll, HPdOl), TPdOU U 12 

EXTERNAL VISCCS U 13 

EXTERNAL PRANDTL U 14 

DO 1 1=1, N U 15 

TPd)=Td) U 16 

CALL OIMPROP (P(I),H(I),PPm,HP(m U 17 

EMUd l = VISCOStTPm,PP(I M U 18 

EMU( t l=EMU(I)/EMUSHOK U 19 

PRANd l=PRANDTL(TP( 11 ,PPd M U 20 

I CONTINUE U 21 

RETURN U 22 

END U 23- 
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FUNCTION RHOSY ( P ,H , AFC t AFN, AFO, TGUES J V 

V 

P ATMOSPHEREStH DEGREES K, RHOSY AMAGATS V 

V 

COMMON /FRSTRM/ H I N , P I N, RHOI NP , U I NP t R 0 V 

V 

INPUTS TO RHOSY P ATMOSPHERES, H DEGREES K V 

INPUTS TO RATRAP P ATMOSPHERES , H CAL/GM V 

OUTPUT FROM RATRAP RHOSY GM/CC V 

V 

CALL RATRAP ( 1 . , 1 , 1 , 1 , TGUES , 0. , P , AFC , AFN, AFO, H/ 1 3.483*A. 19 I , QR, 0 , 2 V 
1, RHOSY, 1,T) V 

TGUES=T V 

RHCSY=RHOSY/RHOINP V 

RETURN V 

END V 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 
16- 


FUNCTION VISCOS (AT,AP) W 

COMMON /HANSEN/ T ABT ( 31 I , T ABPLOG ( 6 ) , T A BZH ( 186 » , T ABZ ( 1 86 ) , T ABPR 1 1 86 W 
1),TEMURAT(1861 W 

EXTERNAL RHOSY W 

INPUT T TEMPERATURE, DEG. K W 

P PRESSURE,ATM H 

OUTPUT VISCOSITY, GM/CM/SEC W 

T=AT W 

P=AP W 

PLCG*AL0G10(Pl W 

CALL DISCOT ( T , PLOG, TA6T ,TEMUR AT , TABPLOG, 21 , 186, 6, V I SCOS » W 

IF (T.LE. 15000. » GO TO 1 W 

WRITE (6,21 T H 

VISC0S=VISC0S*(l.462E-5)*SQRT(TI/(l.*ll2./TJ H 

RETURN W 

W 

M 

FORMAT (6X,25HTEMPERATURE TOO LARGE. T=E16.6, 5H0EG.K ) W 

END W 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19- 


FUNCTION PRANCTL (AT,AP) X 

COMMON /HANSEN/ T ABT ( 31 ) , TABPLOG ( 6 I , T ABZH( 186 1 , TABZ ( 1 86 ) , TABPR ( 1 86 X 
1) ,TEMURAT( 186) * X 

INPUT T TEMPERATURE, DEG. K X 

P PRESSURE ATMOS X 

T=AT X 

P = AP X 

PLCG = AL0G10(P ) X 

CALL DISCOT (T , PLOG,T ABT , TABPR ,T ABPLOG , 21 , 186, 6,PRAN) X 

PRANDTL=PRAN X 

IF (T.LE. 15000. ) GO TO 1 X 

WRITE (6,2) T X 

RETURN X 

X 

X 

FORMAT (6X,25HTEMPERATURE TOO LARGE. T=E16.6, 5HDEG.K ) X 

END X 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 
17 - 
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BLOCK DATA Y 1 

COKMON /HANSEN/ T ABT ( 3 1 » , T AB PLOG ( 6 » , T A BZH { 1 86 I f T ABZ ( 1 86 I , T ABPR < 1 86 Y 2 

1 ) ,TEMURAT{ 186) Y 3 

COPMON /COINGll/ TSTR(18 I.QMEGAims) Y 4 

T TEMPERATURE DEGREES K Y 5 

P PRESSUREt ATM Y 6 

DATA TABT/100., 500. ,1000. ,1500., 2000. t2500., 3000. ,3500., 4000. t4500 Y 7 

1 . . 5000. . 5500. .6000. .6500. . 7000. • 7500. , 800 ). , 8500. ,9000. , 9500. , 1000 Y 8 

20. . 10500. .11000. .11500.. 12 JOO.,12500. ,13000. ,13500. , 14000. , 14500. , Y 9 

315000./ Y 10 

DATA TABPL0G/-4.,-3.,-2.,-l.,0.,l./ Y 11 


DATA TABZH/3. 52, 3. 52, 3. 65, 3. 80, 4. 41, 8. 02, 8. 19, 8. 03, 9. 82, 16. 80, 23. 4 Y 12 
16,23.58,22.54,21.93,22.29,24.26,28.65,35.75,43.74,49.15,50.96,50.6 Y 13 
24,49.48,48.07,46.64,45.26,43.97,42.76,41.64,40.58,39.60,3.52,3.52, Y 14 

33.65. 3. 80.4.07.6.16.8.02 .7.77.8.09.10.55.16.68.21.58.21 .97.21.24.2 Y 15 

40.69.20.72.21.65.23.85.27.66.33.00. 38.79.43.28.45.57.46.09.45.64.4 Y 16 

54.74,43.69,42.61,41.55,40.53,39.57,2*3.52,3.65,3.80,3.97,4.81,7.13 Y 17 
6,7.62,7.53,8.14,10.48,15. 14,19.30,20.35,20.01,19.54,19.34,19.60,20 Y 18 
7.49,22. 17,24.78,28.28,32.31,36.13,39.01,40.66,41.26,41.17,40.69,40 Y 19 
8.01,39.24,2*3.52,3.65,3.80,3.93,4.27,5.55,7.08,7.28,7.33,7.96,9.73 Y 20 
9,12.93,16.46,18.34,18.66,18.43,18.17,18.09,18.29,18.85,19.84,21.31 Y 21 
$,23. 28, 25.69,28.36, 30. 99 , 3 3. 27, 34.97, 36.02, 36. 53 ,2*3 . 52, 3.65, 3. 80, Y 22 
$3.92,4.09,4.61,5.75,6.74,6.98,7.10,7.58,8.70, 10.64,13.20,15.48,16. Y 23 
$73,17.09, 17.04,16.91,16.84,16.90,17.13,17.57,18.24,19.16,20.32,21. Y 24 
$72,23.29,24.98,26.66,2*3.52,3.65,3.80,3.92,4.03,4.25,4.75,5.56,6.2 Y 25 
$9,6.62,6.80,7.11,7.72,8.76,10.24,11.99,13.63,14.79,15.40,15.61,15. ,V 26 
$64,15.60,15.58,15.62,15.74,15.96,16.28,16.71,17.26,17.92/ Y 27 

DATA TABZ/4*!., 1.016, 1.163, 1.2, 1.211, 1.287, 1.577, 1.91, 1.99, 2. 008, 2 Y 28 

1.032.2.088.2.21.2.446.2.826.3.282.3.645.3.843.3.932.3.969.3.985.3. Y 29 
2993,3.996,3.998,2*3.999,2*4. ,4*1., 1.005,1.088, 1.192,1.203,1.228,1. Y 30 

3337.1.622.1.898. 1.983.2.006.2.027.2.067.2.144.2.284.2.51.2.832.3.2 Y 31 
402,3.526,3.745,3.867,3.931,3.963,3.979,3.988,3.993, 3.996,3.997,4*1 Y 32 

5. . 1.002.1.033.1.149.1.197.1.208.1.245.1.359.1.599.1.849.1.961.1.99 Y 33 

67.2.017.2.044.2.09.2.166.2.286.2.462.2.700.2.983.3.272.3.52.3.7.3. Y 34 

781 8, 3. 889, 3. 932, 3. 957, 3. 973, 4*1., 1.00 1,1. 01 1,1. 072, 1.167, 1.198, 1.2 Y 35 
813,1.252,1.348, 1.529,1.752,1.904,1.971,2.001,2.023,2.05,2.09,2.149 Y 36 
9,2.234,2.351,2.505,2.694,2.91,3.135,3.347,3.527,3.667,3.769,5*1. ,1 Y 37 
$.004,1.026,1.092,1.165,1.196, 1.214,1.248,1.316,1.437,1.607,1.778,1 Y 38 
$.896,1.959,1.993,2.018,2.042,2.071,2.111,2.163,2.232,2.318,2.426,2 Y 39 
$.553,2.7,2.861, 3. 028, 5*1 ., 1. 001, 1.009 , 1.035 , 1.089, 1 . 149, 1 . 186 , 1. 20 Y 40 
$8,1.235,1.279,1.351,1.457,1.59,1.727,1.838,1.914,1.962,1.993,2.018 Y 41 
$,2.042,2.067,2.098,2.135,2.18,2.233,2.297,2.372/ Y 42 

DATA TABPR/2*. 738, .756, .767, .614, .771, .714, .606, .587, .764, .993, .87 Y 43 

11. . 4 55, .392, .361, .342, . 322 , . 279, .2 , .1 14 , . 0576, .0314, .021 3,. 0167, .0 Y 44 

2143. . 0129. .0121. .011.. 0108.. 01 09..01 1.2*. 738,. 756,. 767, .668, .654,. Y 45 

3745. . 658. . 58..611 . . 799 . . 989 . . 891 . .464 . .404. . 371 . .35 1 . . 335 . . 316 . . 27 Y 46 

49.. 2 16. .14 5.. 0877. .0524. .0346 . .0238 .. 019. .0162 . .0149 .. 013. . 012 . 2*. Y 47 

5738. . 75 6.. 767.. 7 24.. 611.. 7 4.. 737.. 619.. 578. .624.. 78 5.. 969.. 9 55. .83 Y 48 

6. . 424. .387. .363. .348. .336. .319. .295. .254. .201. .146. .101. .0688. .047 Y 49 

7.. 03 45. .0245. .0129. 2*. 73 8, .756,. 767,. 766, .645, .636, .744, . 759 , .61 , . Y 50 

8 581. . 617. .73 6. .906. .986. . 969. . 648 ..41 1 . . 332 . . 364. .348. . 339. . 327 . .3 Y 51 

9 12.. 292.. 263.. 227, .185,. 144,. 0986,. 08 19, 2*. 738,. 756,. 767,. 773,. 696 Y 52 

$,.627,. 660,. 76 2,. 7 52,. 611,. 58 3,. 602,. 6 73,. 79 6,. 927,. 983,. 943,. 807, Y 53 
$.497, .429, .404, .382, .369, .355, .343, .333, .319, .302, .277, .253, 2*. 738 Y 54 
$,.756,. 767,. 773,. 751,. 680, .6 31,. 662,. 743,. 767,. 62, 2*. 592,. 62,. 688, Y 55 
$.788, .891, .961, .966, .872, .532,. 463,. 434,. 412,. 396, .383, .369, .36, .3 Y 56 
$49, .341/ . Y 57 

DATA TEMURAT/7*!., 1.011, 1.032, 1.096, 1.181, 1.227, 1.256,1. 271, 1.264, Y 58 

11. 2 1.1. 072.. 826. .517.. 261.. 118.. 055.. 029.. 01 8.. 012.. 009.. 008. 2*. 00 Y 59 
27, 2*. 00 8, 7*1., 1.01,1.024,1.055,1.128,1.209,1.257,1.286,1.303,1.307 Y 60 

3.1.28.1.207. 1.068.. 853.. 59 5.. 361.. 2.. 108.. 063.. 036.. 024.. 01 8.. 01 5, Y 61 

4. 013.. 012. 7*1. ,1.01,1.022,1.038,1.074,1.146,1.228,1.276,1.317,1.33 Y 62 

57. 1.347. 1.343. 1.314. 1.2 5 1.1. 143.. 9 83.. 7 82.. 571.. 387.. 249.. 158..!,. Y 63 

6067. . 042. .016.7*1. , 1.006,1.02,1.033,1.051,1.086,1.148,1.229, 1.294, Y 64 

71.332.1.371. 1.386. 1.396. 1.39 3. 1.375. 1.335. 1.267. 1.168. 1.04.. 881. .7 Y 65 

8 11.. 547.. 408.. 262.. 212. 7*1., 1.003, 1.016, 1.029, 1.043, 1.06, 1.09, 1.13 Y 66 

99,1.208,1.283,1.342,1.386,1.425,1.438,1.445,1.448,1.442,1.424,1.39 Y 67 
$4, 1. 342, 1.274, 1.187,1.082, .94,. 828, 7*1. ,1.001,1.008,1.022,1.036,1. Y 68 
$052, 1.067, 1.09, 1.124,1.175,1.238,1.307,1.368,1.418,1.468,1.496,1.5 Y 69 
$01,1.511,1.52,1.516,1.508,1.492,1.468,1.415,1.387/ Y 70 
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DATA TSTR/5. ,6. ,7.f 8.,9. ,10. t 20. t30.,<iQ.,5C.t60.,70., 80. t90. .100., 7 71 

12)0. ,3-)C. ,400./ Y 72 

DATA OMEGAll/. 8422, .8124, .7896, .7712, .7556, .7424, .6640, .6232, 0.596 Y 73 
10,0.5756,0.5596,0.5464,0.5352,0.5256,0.5170,0.4644,0.4360,0.4170/ Y 74 
END Y 75- 


SUBRQUTINE Ol^'PROP ( AP , AH, PPR I HE ,HPR I ME I 2 1 

COMMON /FRSTRM/ H I N , P 1 N , RHO I NP ,U I NP , R fi I 2 

L 3 

INPUTS TO DIMPRQP ARE DIMENSIONLESS I 4 

OUTPUTS P ATOMOSPHERES ,H ERGS/GN I 5 

I 6 

PPRIME=AP*RHOlNP*UINP*UINP/( 1.01325E+Q6I I 7 

HPRIME=AH*UINP*UTNP ' I 8 

RETURN 2 9 

ENC 2 10- 


SUBROUTINE TINT (N,DX,A,8) 

AA 

1 

THIS SUBROUTINE INTEGRATES A TABLE TO GIVE A TABLE 

AA 

2 

DIMENSION A( 1021 , 6( 102t 

AA 

3 

NHAF=N/2 

AA 

4 

IEVEN=0 

AA 

5 

IF (2«NHAF.EQ.NI IEVEN=1 

AA 

6 

IF (lEVEN.EO.l) NHAF=NHAF-1 

AA 

7 

B(1I=0. 

AA 

8 

DO 1 I=1,NHAF 

AA 

9 

J = 2*I 

AA 

10 

B(JI = B( J-1 M-.08333333*DX*I 5.*AL J-l)*8 .♦AI J»-A( J+1 > » 

AA 

11 

J=J+1 

AA 

12 

B( J) = B( J-2» + .33333333*DX*( A( J-2) *4.*A ( J-1I*A( J) » 

AA 

13 

CONTINUE 

AA 

14 

IF ( lEVEN.EQ.O) GO TO 2 

AA 

15 

J=J*1 

AA 

16 

B( JI = B( J-1)*.08333333*DX*(-A(J-2I*8.*AJJ-1)*5.*AU1 ) 

AA 

17 

RETURN 

AA 

18 

END 

AA 

19- 
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SUBROUTINE GET REDV 

DIMENSION PALM ( 100 ». P ALC ( IDO ItPALNI 1001. PALO (100) 


DIMENSION 

A(6,6) 

,A12(16,7) 

,A22(7,7) 

, ALP(6) 

,ALPT( 16,5) 

.FEMP0040 

IBMTC 16) 

,C(6,16) 

,CH(16,2) 

,CP(16) 

,DEGI(6) 

,H(16) 

,FEMP0050 

2JAT116,5) 

,JPH(16) 

,KAT(6) 

,K00E( 16) 

,RA(16,2) 

,R6(16,2) 

,FEMP0060 

3RC(16,2) 

,R0{16,2) 

,RD1( 16,2) 

,RE(16,2) 

,RE1(16,2) 

,SD( 16) 

,FEHP0070 

4TB(3) 

,TU(16,2) 

,TU2( 16,2) 

,VN(17) 

,VNE(16) 

,VNT(17) 

.FEMP0080 

5VNUC 16,6) 

,W26(6) 

,W3(16) 

,Y(16) 

,RF(16,2) 

,RC1(16,2) 

FEMP0090 

DIMENSION 

SAC (200), 

SAH(200) 





COMMON /LOKHEED/ 






A A 

,A12 

,A22 ,ALP 

,ALPT 

,BHT ,C 

, CH ,- 

FEMPOlOO 

ICP ,OEGI ,H 

,HS ,IG 

,IGMS 

,IGMSP ,IGP 

,ION , 

FEMPOllO 

2IS ,ISP 

, ISPNG 

, ISPNGP, ISPNG2, JAT 

,JPH ,KAT 

,KOOE , 

FEMP0120 

3N ,NG 

,NP 

.PRESS ,RA 

,RB 

,RC ,RD 

,R01 , 

FEMP0130 

4RE ,RE1 

,S0 

,TU ,TU2 

,VN 

,VNE ,VNT 

, VNU , 

FEMP0140 

5W26 ,W27 

,W3 

,Y ,RHO 

,WM 

,SYU ,RF 

.RCl.TB 

FEMP0150 


COMMON/RA0/YY(100),TEE(100l,FHV(20) ,NHViNT,C2, lY 
C0MMON/RAD/QRYPC(100l ,QRYPL( 1001 
COMMON/ RAD/ XNN (14,100). XMOL 

COMMON/ RAD/ N1 HVC , FHVC ( 50 ) , AHV( 50 ) , AHWL ( 23 ) 

COMMON /RAD/C1,C3,C4,FLG,C5,FLG1 
COMMON /RAD/ YDELT , DELTA, FLl , FL2 

COMMON/ RAO/GEE(8) , EPS ( 8 ) , NU( 20 ) . NO ( 70 ) , FF ( 70 ) .GAMP (70), 
1W0L(20),FHVM(20),FHVP(20) 

COMMQN/CIONCL/F(100,10),F2(100,10),HVL(70),EP,K2,K1, IFL, lYCON, 
lIOl.WMI , BIJ( 100. 10 ) ,GMIN( 100.10) ,GPLU( 100, 10) ,IAED 
DIMENSION NICN(IOO) 

DIMENSION PRES(IOO) 

DIMENSION HH(IOO) 

DIMENSION 2BLK(1157) 

DIMENSION QRYP(IOO) 

EQUIVALENCE (ZBLK.A) 

XM0L=1. 

00 1157 1=1,1157 
1157 Z8LK(n=0. 

NHV=18 

00 998 1=1,17 
998 VNT(I)=0. 

C4=1.273 

G=l. 

100 F0RMAT(6E12.1 ) 

INPUT TABLES I, I I 

REA0(5,100) ( GEE( I ) . 1 = 1, 8) 

READ(5,100) (EPS( I ) ,1=1,8) 

READ(5,100) (FHVM( I ),I = 1,NHV) 

READ(5,100) (FHVPU ),I = 1,NHV) 

READ(5,100) ( FHVC I ) ,I = 1,NHV) 

READ(5,100) ( WOL( I ) , 1=1 ,NHV) 

READ(5,101) C NU( I ),I = 1,NHV) 

101 F0RMATC18I2) 

IS=0 

DO 1 1=1, NHV 
1 IS=1S+NU(I) 

READ! 5, 102) (NOC I ),HVL( I ),FF( I ) ,GAMP( 1 1 ,1 = 1, IS) 

102 FORMAT! II, 11X,3E12.1) 

INPUT HV VALUES AND AK FACTORS FOR CGNTM CALCULATION 

READ(5,115) NIHVC 
115 F0RMATC24I3) 

READCS.lOO) (FHVC( I) ,I = 1,NIHVC) 

READ(5,100) ( AHVC I ) , I=1,NIHVC) 

READ (5,100) (AMVLC I ),I = 1,NHV) 

RETURN 

END 
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SUOPOUTIM?: WATP&P ( OALTA» NIC, NICN, NX. TFX, YX, PRFS* PALC, 
1PALN,OALO,HH,QPYP. IPAD, IOPT,PO,KI)T,TEFE> 

OIMFNSION PALHdO M .PALCdOO) ,PALNdOO> .PALOdOO) 


DIMENSION 

A ( 6 , 6 ) 

,A12(.16,7) 

,A22(7,7) ,ALP(6) 

,ALPT(16,5) 

,FEMPO04O 

IPmt ( 16) 

,C (6, 16) 

,CM(16,2) 

,CP(16) ,nEGI(6) 

,H(16) 

.FEMPOOSO 

2JAT (16,S) 

, JPM (16) 

,KAT (6) 

,K0DEd6) ,RAd6.2) 

,RB (16,2) 

,FEMP0060 

3RC(1G,2) 

,R0 ( 16,^) 

,RD1 (16,2) 

,RE(16,2) ,RE1(16,2) 

,SD(16) 

.FEMP0070 

4T0 (3) 

, TIJ(16,^) 

,TU2(16,2) 

,VN(17) ,VNE(16) 

,VNT(17) 

.FEMPOOSO 

SVNU(10,6) 

,W26 (6) 

,\«)3(16) 

,Y(16) ,RF(16,2) 

,RC1 (16,2) 

FEMP0090 

dimension 

SAC (200) . 

SAH(200) 





common /lokheeo/ 






A A 

,A12 

,A2? ,ALP 

,alpt ,RMT 

,C 

,CH , 

FEMPOlOO 

ICP ,DEGI ,H 

,HS ,IG 

,IGMS .IGPSP 

, IGP 

,ION , 

FEMPOllO 

2TS ,1SP 

,1SP 'G 

,ISPNGP,IGPNG2,JAT ,JPH 

,KAT 

.KODE , 

FEMP0120 

3n ,ng 

,NP 

.PRESS ,RA 

,PR ,PC 

.RO 

,R01 , 

FEMP0130 

4PE ,RE1 

,S0 

,TU ,TU2 

,VN ,VNF 

, VMT 

,VNU , 

FEMP0140 

5W26 ,W27 

,W3 

,Y ,RHO 

,WM ,SYU 

,RF 

.RCl.TB 

FEMP0150 


COMMON/PAO/YY dooi .TEEdOO) ,EHV<P0> , NHV , NY , C? , I Y 
COmmON/RAO/ORYPC ( I nO) .QRYPL dOO) 

COYMOM/RAH/XNNd A. ) om , XMOL 

CO«MOM/RAD/ NIHVC-EHVC(SO) ,AHV(50) ,AHVL(20» 

COMMON /PAn/CI ,C3,CA,ELG,CS,ELr,l 
COMMON /RAD/ YQELT, delta, ELI ,EL2 

COMMON/RAD/GE'^ ) .PPS(H) ,NII(20) ,N0(70) ,FF(70) ,GAMP(70> . 
1wol(2'’1 ,FHVm(20) hVP(20) 

COMMON/C toncl/f d •■'0, 1 0) ,f 2 don, 1 0) ,hvl (70» ,ep,k2,ki . tfl, i ycon, 
ITOl ,WMT,RIJ (100,1 M ,GMIN (100,10) , GPLU (1 00 , 1 0 ) ,IAEO 
DTMFNSION NICNd 0) 

0Tm='N=)TOKi PRESdO'M 
OTmFMSION HH(IOO) 

OImfaj<;ION /HLOdl-’M 
dimension ORYPdO. ) 
dimension YXdOO) 
dimension POdOOl 
DJMFNSION Kl?®!) 

DTmFNSION XNC2H('.0) ,XNC3H(A0) ,XNCAH(40) .XNHCA'(AO) ,XNC?H2(A0) , 
IXNHP(AO) 

DTMFNSION TEEE(IO') 

POUT valence (7PLK.A) 
delta=oai;ta 

NY = ^!X 

00 ?00 1=1, NY 
200 YY(T)=YX(I) 

FL1= 2. 

EL2= 2. 

ELG= 0. 

ELG1=1 . 

T6d )=TEX 
L0NGRTT=0 
PUTl 2=100. /12.. 

PUT14=100./14. 

PUT1G=100./1G. 

DO P2'S 1 = 1, NY 
SAC(I)= PALCd) 

SAP(T)=1 .0-(PALC( I ) *PALN ( I ) +OAL0 ( I ) ) 
palC(T)=pUT12*PALC(I) 

DALN( I) =PUT14*PAL n( I ) 

PAL0(n=oUT16«PAL0(n 

PALM ( I) =1 00.-(12.odalC( I) ♦14.»PALN(I) ♦le.oPALOdn 
IF (PALCd ) .LT. 1 .OF-nS) dALH(I)=0.0 
22S CONTINUE 

4 DO 20 1=1, NY 
ppFSS=PRFS(I 1 
IE ( IOPT.EO.2) go to 7 

given t, get H 

T0(2)=Tfld) 

GO TO H 
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APPENDIX A - Continued 


C GTVFN' H, GFT T 

7 Tt3(?>=T8(n*l.? 

HS=HH(T) 

6 Tl=n 
fi’ CONTTMUF 

TF (T8 (?) .LT.GOOn, 1 GO TO 9 

HIGH TFMDF94TURF r^LCULariON 

4L°(?>=PALH(D 
4l_P(G)=P4LC(I) 

4LD (4) =P4LN ( I ) 

4LP(5)=PaL0(I> 

4LD ( 1 ) =ALP (3) ♦ 4lp(4) ♦ ALP(5) «ALP(2) 

IT = ) 

TTXV=-l 

CALL HTFST (ITX,1TXV,1) 

TTXV=ITX 

CALL FFMD (iHFLPffOPT) 

CALL HTFST (ITX.irxv,?) 

TZw = ? 

GO TO in 

LOW tfmpfraturf . 

9 IT5V=-1 
IT=1 

SLD(n=P4LH(t) 

AL°(?)=P4LC(n 
4LD(3)=PALN(I) 

4LD(4)=P4L0(I) 

CALL LTFST (TTSiIT.SV,!) 

ITSV=ITS 

CALL FFMD (IHFLP»(0PT) 

CALL LTFST (ITSflTPv,?) 

T7W=1 

iO IF (lOPT.EO.?) GO TO 32 
T0(?)=TFF(I) 

T8(1)=TR(2) 

GO TO 33 

32 TFF(I)=TR(2) 

M5=HM ( I ) 

Tan ) =TR (2) /I .2 

33 IF (IhPLP.FQ.O) GO TO U 

^RPOP IN FEMP 

WRITE (6t nn) I »PRFSS,HS»TR ( ] ) ,TR<?) 
lio PORM4T ( 1 OHIFEMP Rl.FW » II 0 . 4F20 . 6 ) 

STOP 

RACK TO (12) JUST one TIMF...LOOP COULD OCCIIPF OTHPRWISF 

U TF( I) .NF.O) GO TO 13 
11=1 

GO TO 1? 

CHPCK FOR high OR LOW TEMP CALCULATION 
...THPN COMPUTE N". DENSITIES 

13 IF ( T 7W.PO.?) GO T ' IS 
*NHP(T)=n.n 

XNCPH ( I) =FI.UT (RMT (H) ,W3(0) ,RHO) 

XNC3H ( 1 ) =FLUT (HMT ( 1 0) «W3 ( 1 0) .RHO) 

XNC4H ( I ) =FLUT (BMT <1 1) ,W3 ( 1 1 ) .RHO) 

XNhCM ( I ) =FLUT (RMTM?) .WT(1?) ,RHO) 

XNC?H? ( 1 ) =FLUT (riMT ( 1 3) .W3 ( 1 3) ,RHO) 

XNM ( 1 . I ) =0, 0 

XNN ( ?, I ) = FLUT (H T (6) .W3 (G) .PHO) 
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■ XNN(?*n=0.0 

XNNUi., I ) = rUHT (B- T (P) ,1^1 («) ,RHO> 
xMV(s«n='i.o 

XNM(B»I)= KLUT(H iT(3) .W3(3) ,PH01 
XMM(7. n = FLMT (H'^T (4) ,W3 (4) .RHO) 

XNM (B. n = FL'IT (8 !T (1 S) « W3 ( IS) ,RHO) 

X\)l^(0, r ) = FL'IT (8"T (14) , W3 (14) ,RHO) 

XMM n 0, T) = FUJI ( (?) • w3 (2) »RHn) 

XMvn 1 .1) = FUJI ( 'MT ( 1 ) .W3(l) .RHO) 

XN'M(1?,I)= FLUT( -"T(7) .W3(7) .RHO) 

XHM(1 ■).I)=C.O 

XNIH ( 14. T ) rFLIJT (HMT (S) .W3 (S) .RHO) 
so TO 19 

IS XN\i(l.I)= RLIIT (8 'T (6) . W3 (6) .RHO) 

XNHP( T ) =FLUT (8MT CO .W3 (3) .RHO) 

XNC 2 H ( I ) =o.n 

XNC3H(I)=0.n 

XNC4H(I)=0.0 

xhhCM(I)=0.0 

XNJC?H?{T)=n.n 

XNM(?.I)= FLUT(8 iT(p) ,W3(8) .RHO) 

XMN(3.T)= FLMT ( -imT (7) .W3(7) .RHO) 

XnN(4.T)= FLUT(H-'T(5) .w3(5) .RHO) 

XNM (S. I ) = FUJI (PMT (S) .W3 (6) .RHO) +FLUT (RMT (7) .W3(7) .RHO) ♦ 

IFLUT (8HT ( 1 1 ) .VJ3 ( 1 ) ) .RHO) .XNHP (I) 

Xf\|M (6. I ) =0.0 

XNN(7.T)= PLI)T(ri 'T(9) ,W3(9) .RHO) 

XM\'(P. I ) = FlUT (8 iT ( 1 0) .W3 ( 10) .RHO) 

XHN(9.n=0.0 

XN\O)0.T)=n.n 

XNM ( 11 . T )'= FLUT (^mT(4).W3(41 .RHO) 

XMN(1?.I)'= FLUT( -i>n(?) .W3(2) ,RHO) 

XMM( 1 3. n = FLUT(8MT( 1 I ) .W3 ( n ) .RHO) 

XNN(14,I)=FLIJT(BHT(12) ,W3(12) .RHO) 

SFT NFW OUFSS FOR temp. . .OPPLIFS ONLY IF IOPT=2 

19 TB(?)=TFF(T) 

hh(T)=hS 

RO(T)=RHO 

IF ( KUT.FO.l ) GO TO 20 
PRINT FFMP CALCUL'iTIONS 

WRTTF (6.109) 

109 FORMAT(14H0 PATH LENGTH. 4X. 

ISHPPFSSURF.GX. 1 1HTFMPFRATURE.4X.8HFNTHALPY.SX.7HDFNSITY) 
o'RTTF (6. 1 1 1 ) YY ( I ) .PPFS(T) .TEE ( I) .HS.RHO 
111 FORMATdOEl.l.S) 

XT = 0. 

DO 2 J=1.14 
? XT= XNN(J.n*XT 

XT = XT*XNC2H( 1 ) *XNr.3H ( 1 ) ♦XNC4H ( 1 1 .XNHCNd 1 ♦XNC2H2( 1 ) ♦XNMP (1 ) 

“00 3 J=1.14 
3 X(J)= XNN(J.I)/XT 
X ( IS) =XNO?H ( I) /XT 
X ( 1 6) =XNC3H ( I ) /XT 
X(17)=XNC4H(I)/XT 
X ( 1 B) =XNHCN ( I ) /XT 
X(19)=XNC2h 2(1)/XI 
X(20)=XNHPm/XT 
WRTTF (6.900) ( X ( I) . J=1 . 20) 

900 FORMAT (SX14HM0LF FRACTIONS /6X . ?HO» . 1 OX . IHO . 1 3X . 2HN* . 1 1 X . 1 HN . 1 1 X . 

1 1hf,irx.2H02.12X.phN2 / 6X.2HCO.10X.2HH2.12X.2HC2.10X.2HCN.11X. 

2 1HC.12X.2HC*.12X,bh H/6 X . 3HC2H.9X . 3HC3H. I 2X . 3HC4H. 1 0 X , 3HHCN.9X . 

3 4hC?H2.9X.2HH*/7M 3.5/7E13.S/6E13.5) 

20 COMTI-JUF 

DO 1 T=1.NV 
1 TEEE ( T ) =TFF ( I ) 
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APPENDIX A - Concluded 


IF (T‘?an.F.Q.n> go to ip 3R 
TN'PIIT 4K F4CT0RS -OR LINr CALCULATIONS 
MW=M4X0 (NHV*NIHVC) 

IPS FORMAT (n P.3 , f 16. n 

1 ?7 F0RM4T(34X,F1?.3,- 16.31 
V5R F0RMAT(2(FI2.3,Fl-^.3n 
IF (FLGI .EO. 0. ) GO TO 60 

MULT, factors HY for FLUX CALCULATIONS 

no 61 1=1, NHV 
61 AHVL ( n =AHVL( n *3. 1416 
00 6? 1=1 ,NIHVC 
6? AHV ( n =AHV ( 1 1 *3. I'-l 6 
60 C2=0FLTA«FL1 
Cl =0RLTA«FL? 

C3=C1/3.1416 
C5 = C1 
FP=.nil 

OUTPUT TABLE II 

ipo rORMAT ( !H1 ,60X,6 -itabLF III 
1?1 format (6H1GROUP,-<X,2hHV.12X,3HHV*,11X,3HHV-,10X.1HN,PX. 
13-iwOL , 7X,4HK ( 1 1 , ,5HHV ( 1 1 , 10X,4HF ( 1 1 ,1 1 X,6HGAM (111 
IP? F0RMAT(I4,Frp.3,?f 1 4 . 3 , II ? , F 1 2 . 3 , 1 9 ,F 1 4 . 3 , 1 R?E 1 6 . 2 ) 

IF (LONGRIT.FO.O) go TO 1066 
WRITE (6,1211 

1 ?4 RORMAT (6BX, I9,F14. 3,4X,1DF12.2, 1PF16.2) 

TC1=0 

00 30 1 = 1 »NH\/ 

IC? = ICl ♦ NU(t) 

IC1= TC1*1 

VRITF (6,122) I,F-i\mI), FHVP(I), FHVM(I), NU(I), WOL(I), 
INO (icn ,HVL ( I CD ,FF ( I cn , GAMP (ICl 1 
IF ( TCI .FO. IC2) GO TO 30 
IC3=IC1+1 
00 2S J=IC3,IC2 

pc WRITE (6,124) NO (,I),HVL(J)? FF(J), GAMP(J) 

30 IC1=IC2 
1066 CONTINUE 

00 TOO L=1,NIC 
IY=MICN(L) 

I YC0N=IV 
VOFLT=YY(IY) 

CONTINUM calculation 
C4LL COMTM 

LINE CALCULATION 
CALL LINF 

0RYP(L)=0RYRC(IY) ,0WYPL(TY) 

30 c continue 

IF (F|_G1 .FO.O.O) -RTUPN 
no 434 1=1, NHV 

434 AmVL ( n =AHVL ( I ) /3. 141 6 
00 43S I=1,NIHVC 

435 AHV ( n =AHV ( I ) /3. 1--1 6 
1333 CONTTNUF 

00 AS 1=1 , NY 
DALC,(I)=.12»PALC( M 
oalN(T1=.14*PALN( M 
PALO ( I ) =. 16»PAL0 ( I ) 

6s continue 
return 
e’no 

» = U3 1 
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SAMPLE INPUT 

The input for a sample case is shown in this appendix as follows; 


$NAH4 



$NAM5 


- 

N 

= 

21. 

IRFLUX 

a 

1. 

lOlM 

= 

2. 

lUPOATE 

m 

If 

IRODAMP 

= 

0. 

ALPINJC 

ae 

a 

o 

• 

o 

HDAMP 

= 

0«7£«^00* 

ALPINJO 

B 

0 . 22 E«^ 06 t 

ROOAMP 

= 

0.7E-»00t 

ALPINJH 

a 

o 

• 

o 

£P 

= 

0. lE-Ql, 

ALPINJN 

a 

0.78E+00. 

ERGV 

= 

O.lE-01, 

020 P2C 

a 

0.25E^01* 

EA 

= 

O.lE-01. 

AS 

a 

0.122E«01» 

EPSX 

= 

O.lE-01, 

PINP 

a 

0.l97E403f 

ERO 

* 

0.2E-01, 

RHOINP 

M 

0.272e>06» 

MAXTIME 

= 

2. 

UINP 

m 

0.1525E+C7, 

NIC 

= 

11. 

RB 

m 

0.3048E«C3f 

CURV 

= 

O.lE^Ol. 

ROVW 

a 

-O.lE+00, 

CURVZ 

= 

0« 12E+01, 
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a 

0.28E>01t 

BS 


0. 1E>01. 

ALPINFN 
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0.78E-*^00» 
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= 

0.2E+01, 

ALPINFO 
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0.22E^O0» 
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= 0.195E + n3. 0.183E*03» 0.17E + 03. 0.157E + 03r 0 . 1 43E + 03 ,J 3E^,03 » , 

0.U2E+03, 0.93E+02. 0.72E+02. 0.45E+02* Q.32E+02t 0.29E*02* 

0.27E+02* 0.271E*02, 0.257E+02t 0.245E+02, 0.234E+02, 

■ 0.223E + 02* 0.213E*02. 0.194E + 02, 0. 1704925E*02» . I » .1 . I , l.v I» I» 

I » I » If If If If If If If If If If If If If If If If If If If If If 

If If If If If If If If If If If If If If If If If If I. If If If If 

If If If If If IfJf If If If If If If If If If If If If If If If If 
Iflflflflf 

90V = -O.lE+OOf -0.99E-0If -0.976E-01f -0.95E-0If -0.909E-01f -0.852E-0If 

-0.777E-01f -0,68E-01f -0.548E-01f -0.39IE-01f -0.I49E-01f 
0.294E-01f 0.912E-01f O.I65E*OOf 0.252Ef00f O.SSE+OOf 0.46E+00f 

0.58E*00f 0.71E*OOf 0.851E + 00f O.lEf-Olf If If If If If If If If 

If If If If If I f I f If If If If If If If If If I . If I f I f I f If If 

If If If I« If If If If If If If If If If If If If If If If If If If 

Iflf If If If If If If If If If If If If If If If If If If If If If 
If If If 

■fiEND 
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_ A = O. Of 0.9E-02f 0.2E-01f 0.33E-01f 0.45E-01f 0.55E-01f 0.73E-01f 

- 0/9E-0rf ” d/i iTE*O0f ■ d. 16E*00f O.29E + 00f U.4E + 00f O.48E*00f 
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_P =_ 0.976E+O0f 0.97t.E*0Uf 0.976E*00f 0.976E+00f 0.976E*00f 

d.976Ef00f 0.976E + O0f 0.976Ef-OOf 0.976Ef00f 0.97bE*Of)f 

0.976E+dOf 0.97bE*00f 0.9?6EfOOf 0.976Ef00f 0.975E+00f 
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^ If If If If If If If If If If If If If If If If If If If If I* If If 

. Iflf If If If If If Iflf If If If If If If If If If If If If If If 
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SEND 
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APPENDIX C 


SAMPLE OUTPUT 

In this appendix, the output for the calculation described in appendix B is given. 
First, identification of the case is printed out. After completion of the iterations the 
flow profiles are given, after which the equilibrium composition profiles are output by 
RATRAP. Only a sample of the composition profiles are shown here. Finally, the pro- 
file of the radiation flux vector is printed out. The output is given as follows; 
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DELTA- 1 
ETAU ) 

0. 

S.OOOOOOE-02 
l.OOOOOOE-01 
1.500000E-01 
2.000000E-01 
2. SOUOOOE-Ol 
3.0Q0000E-U1 
3.500000E-01 
A.OOOOOOE-01 
A.500000E-01 
5.0U0000E-01 
S.SOOOOOE-Ol 
6.000000E-01 
6.500000E-01 
7.000000E-01 
7.500UOOE-01 
S.OOOODOE-Ol 
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7.500000E-01 
8.00000UE-01 
8.500Ua0E-01 
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•43321101E^00 
Y( I ) 

0 . 

3.82I030E-04 
7*962S50E-O4 
l*241459E-03 
1.723479E-03 
2.250103E-03 
2,642980E'03 
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5«621403E-03 
7,506546E-03 
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l*235369E-02 
I.500765E-02 
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3«3B3703E>02 
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l.OOOOOOE^OO 
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1, 2764976*04 
1,3007046*04 
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A(I I 
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9.5362 056-03 
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7,2305206-02 
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6,2335646-04 
6,3402086-04 
6.2570306-04 
5.9760216-04 
5,3743476-04 
2,3341056-04 
1,9264096-03 
4,3065686-03 
7.3603736-03 
1,1112216-02 
1,5636436-02 
2,0988536-02 
2,7275956-02 
3.4617596-02 
4,3367006-02 
5,8603776-02 
R0( 1 1 

l.9515276*‘02 
1,8044906*02 
1,6627456*02 
1,5489806*02 
1,4113306*02 
1,2860556*02 
1, 1191896*02 
9,46406 lE^Ol 
7.437146E^01 
4,6381996*01 
3, 1780056*01 
2,9573436*01 
2,775983£*01 
2,6294876*01 
2,5075366*01 
2,3943266*01 
2,2889666*01 
2.1859356*01 
2,0840016*01 
1. 9674996*01 
1,7059906*01 
AMFOl I } 

2,2000006-01 
2.2C000C6-01 

2,2000006-01 

2,2000006-01 

2,2000006-01 
2,2C00CC6-01 

2,2000006-01 

2,2000006-01 

2,2000006-01 
2.2C0000E-01 

2,2000006-01 

2,2000006-01 
2*2000V0t-Ql 
2,20000C6-01 

2.2000006-01 

2,2000006-01 

2, 2000006-01 

2,2000006-01 

2,2000006-01 

2,2000006-01 

2,2000006-01 


04 

6NUII) 

6,4369176-01 
6,6280986-01 
6,8359266-01 
7, 1267516-01 
7.4941846-01 
7.9116036-01 
8,4823586-01 
9*2312936-01 
1.0404246*00 
1,3663116*00 
1,5984286*00 
1,5823746*00 
1,5501406*00 
1,5123316*00 
1,4734446*00 
1,4315496*00 
1,3864376*00 
1,3391136*00 
1,2909056*00 
1,2454666*00 
9.99365BE-01 
XNRQU) 
1.9514126*02 
1 , 8112626*02 
1,6666266*02 
1,5401676*02 
1,3946646*02 
1,2699166*02 
1,1085156*02 
9,3899616*01 
7,390891£*01 
4.617C98E*01 
3,1826106*01 
2,964K36*0l 
2,7787876*01 
2.6297576*01 
2,5065926*01 
2, 3928616*01 
2. 2874556*01 
2.1836536*01 
2,0808856*01 
1,9522436*01 
1,7061986*01 
AMFHU ) 
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0 , 

0. 

0 , 

0. 
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PAANin 
7. 0790646-01 
7,4091876-01 
7.6025146-01 
7.5100836-01 
6,9025366-01 

6, 0595316-01 
5,9009886-01 
6,3467866-01 
7,6300116-01 
6,4896536-01 
3,6195226-01 
3,4750656-01 
3.3563646-01 
3,2628376-01 
3,1676966-01 
3^1070116-01 
3,0184666-01 
2. 9275296-01 
2,8365956-01 
2. 7523806-01 
2, 3066486-01 

HU ) 

2, 8000006-02 
3, 1381366-02 
3,45811 IE-02 

3, 7705676-02 
4, 144261E-02 
4.6599786-02 
5.617058E-02 

7. 372470E-02 
I.O86976E-01 
1, 659130E-01 
2.5478156-01 
2.757263E-01 
2,969440E-01 
3, 1655186-01 
3,3460646-01 
3. 5279416-01 
3, 7092596-01 
3. 8998586-01 
4,0994606-01 
4.3675586-01 
4,9924416-01 

AFORUl 
9. 987 159 6-01 
9. 9790706-01 
9.9645626-01 
9. 9387456-01 
9. 8931426-01 
9. 0123546-01 
9.6661696-01 
9. 3873666-01 
8. 8113356-01 
7.4240726-01 

4. 1831616-01 
1,1713496-01 
1. 8976756-02 
I. 9608806-03 
1. 3885506-04 
7. 0893096-06 
0 . 

0 . 

0 . 

0 . 

0. 


YOYSd ) 

0, 

1,0130166-02 
2,1111066-02 
3.2913036-02 
4.5692166-02 
5. 9653786-02 
7.5371916-02 
9.3625486-02 
1.1616346-01 
1.4903236-01 
1.9906346-01 
2.6112746-01 
3.2751596-01 
3,9787656-01 
4.7187866-01 
5,4938386-01 
6,3049836-01 
7.1536246-01 
6.0428486-01 
8.9707336-01 
1.0000006*00 
XNH< I ) 

2. 8000006-02 
3.0812466-02 
3,4258516-02 
3.8422606-02 
4,2987676-02 
4. 8775076-02 
5.8896546-02 
7.6940936-02 
I. 1207026-01 
1.8796346-01 
2. 5346196-01 
2.7333276-01 
2.9577266-01 
3.1638656-01 
3.3507026-01 
3. 5360226-01 
3.7181996-01 
3.9147966-01 
4. 1201716-01 
4.4891696-01 
4.9937286-01 
AONC( 1 ) 

1. 2640706-03 
2.0929896-03 
3.5417936-03 
6.1255056-03 
1.0665616-02 
1,8764616-02 
3, 3381056-02 
6,1263446-02 
1, 1886656-01 
2.5759286-01 
5, 8168396-01 
8,8286516-01 
9»8l0233B‘0l 
9,9603916-01 
9.9986116-01 
9.9999296-01 
1.0000006*00 
1,0000006*00 
1,0000006*00 
1,0000006*00 
1,0000006*00 



APPENDIX C - Continued 


PATH LENGTH 

PRESSURE 

TEMPERATURE 

ENTHALPY 

DENSITY 



0. 

6. 09388E-01 

3. 61398E+U3 

1.55566E+03 

5.30A59E-05 



MOLE FRACTIONS 






0+ 

0 

N+ 

N 

E 

02 

N2 

CO 

F2 

C2 

CN 

C 

C+ 

H 

C2H 

C3H 

CAH 

HCN 

C2H2 

H + 


0. 

2.07175E-QL 

0. 

3.99353E-0A 

0. 

7.38121E-02 

7.1861AE-01 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 


path length 

PRESSURE 

temperature 

ENTHALPY 

DENSITY 



1.013C2E-02 

6.Q9386E-01 

3.79647E+03 

1.7A352E + 03 

A.92326E-05 



MOLE FRACTIONS 






0+ 

0 

N* 

N 

E 

02 

N2 

CO 

H2 

C2 

CN 

C 

C+ 

H 

C2H 

C3H 

CAH 

HCN 

C2H2 

H + 


0. 

2.5I558E-01 

0. 

8.56953E-0A 

0. 

4.7182AE-02 

7.00AC2E-01 

0. 

0. 

0. 

0. 

C. 

0. 

0. 

0. 

0. 

0. 

U. 

0. 

0. 


PATH LENGTH 

PRESSURE 

TEMPERATURE 

ENTHALPY 

DENSITY 



2.LUUE-02 

t.QS384£-01 

3. 99909E+03 

1.92130E+03 

A.57768E-05 



MOLE fractions 






0+ 

0 

N+ 

N 

E 

02 

N2 

CO 

H2 

C2 

CN 

C 

C* 

H 

C2H 

C3H 

CAH 

HCN 

C2H2 

H + 


0. 

2.e6503£-(U 

0. 

i.851i8E-03 

0. 

2. 6153 3E-02 

6.85AS3E-0I 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 


PATH LENGTH 

PRESSURE 

temperature 

ENTHALPY 

DENSITY 


- 

3. 29130E-02 

6. C9383E-0L 

A.29857E+03 

2.09A9OE+03 

A.185UE-05 



MOLE FRACTIONS 






0+ 

0 

N+ 

N 

E 

02 

N2 

CO 

H2 

C2 

CN 

C 

C + 

H 

C2H 

C3H 

CAH 

HCN 

C2H2 

Ht 


0. 

3. 12849E-01 

0. 

5. 09863E-03 

0. 

1.00510E-02 

6.72001E-01 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 


PATH LENGTH 

PRESSURE 

temperature 

ENTHALPY 

density 



4.56922E-02 

6.09382E-01 

A.68828E+03 

2.30252E+03 

3.78972E-05 



HOLE FRACTIONS 






0+ 

0 

N* 

N 

E 

02 

N2 

CO 

H2 

C2 

CN 

C 

C + 

H 

C2H 

C3H 

CAH 

HCN 

C2H2 

H+ 


0. 

3. 22968E-01 

U. 

1. 58038 E-02 

0. 

2.93056E-03 

6.58298E-01 

U. 

C. 

0. 

0. 

0. 

0. 

0. 

U. 

0. 

0. 

0. 

0. 

0. 
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APPENDIX C — Continued 


TOTAL QMINUS 
TOTAL CPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


TOTAL CMINUS 
TOTAL QPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


TOTAL CMINUS 
TOTAL QPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


TOTAL QMINUS 
TOTAL QPLUS 


= 0 . 

= 20.72S6E<-02 


= 56.802SE-09 
= 2i.0752E*02 


= 29.7649E-07 
= 21.28fc9E+02 


= 16.8581E-05 
= 21.0812E+02 


= H.5100E-03 
= 20.68UE+02 


» 77.6169E+00 
= 22.1093E+02 


= 54.7567E+01 
= 24.8172E+02 


= 11.5258E + 02 
= 27.2646E+02 


= 18.3072E*02 
= 28.2948E+02 


= 26.9051E+02 
= 27.4797E+02 


= 41.3077E+02 

= 0 . 



APPENDIX C - Concluded 


ETA(I) Y(l) 

0 . 0 . 

5.000000E-02 3.821030E-04 

l.OOOOOOE-Cl 7.962958E-04 
1.500000E-01 lo241459E-03 

2.000000E-01 1.723479E-03 

2.500000E-01 2.250103E-03 

3.000000E-C1 2.842980E-03 

3.5OQ000E-0i 3.531493E-03 

4.000000E-01 4.381608E-03 

4.500000E-01 5.621403E-03 

5.000000E-01 7.508546E-03 

5o500000E-Cl 9. 849560E-03 

6.000000E-01 1.235369E-02 

6.500000E-01 1.500765E-02 

7.000000E-01 1.779896E-02 

7.50000<JE-C1 2. C72241E-02 

8.000000E-01 2.378199E-02 

8.500000E-01 2.698301E-02 

9.000000E-01 3,033711E-02 

9,500000E-01 3.383703E-02 

l.OOOOOOE+00 3.771936E-02 


YOYSm QRAOU)tW/Sg CM 

0, -3.882058E+03 

1.013016E-02 -3.909761E+03 
2.111106E-02 -3.936354E+03 
3.2913U3E-02 -3. 960952E + 03 
4.569216E-02 -3.982931E+03 
5.9653 78E-02 -4. 004473E+C3 
7.537i91E-02 -4.024465E+03 
9.362548E-02 -4. 076455E+C3 
1.161634E-01 -4.168434E+03 
1.490323E-01 -4. 294697E+03 
1.990634E-01 -4. 504456E+C3 
2.611274£-0i -4.497930E+03 
3.275159E-01 -4.246319E+03 
3.978765E-01 -3. 947319E+03 
4.718786E-01 -3. 512482E+03 
5.493838E-01 -2.925811E+03 
6.304983E-01 -2. i67645E*03 
7.153624E-01 -1.123870E+03 
8.042848E-01 2. 123313E«^02 

8.970733E-01 3.505130E+03 

1. 000000E4-00 8. 8O4495E4-03 
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APPENDIX D 


LANGLEY LIBRARY SUBROUTINE ITRl 

/ 

Language : FORTRAN 

Purpose : To solve the single equation of the form x = f(x) for one real root by the 
Newton- Raphson iteration method. 

Use : CALL ITRl (X, DELTX, FOFX, El, E2, MAXI, ICODE) 

X 

DELTX 

FOFX 
El 
E2 

MAXI 
ICODE 

ICODE = 1: Maximum iteration exceeded. 

ICODE = 2: Derivative = 0. 

Restrictions : A function subprogram with a single argument x must be written by the 
user to evaluate f(x). The name given to the FOFX subprogram must appear in an 
EXTERNAL statement in the calling program. 

Method : The Newton-Raphson iteration technique (ref. (a) of this subroutine) is 
used where 

^n-i-1 = Qn + (1 - ^ f(xn) 


An initial guess supplied by the user. On a normal return to the calling 
program from ITRl, X contains the root. 

An increment supplied by the user so that ^ is a 

reasonable approximation to the derivative of f(x). 

A function subprogram to evaluate f(x). 

Relative error criterion. 

Absolute error criterion. 

A maximum iteration count supplied by the user. 

An integer supplied by ITRl as an error code. This code should be 
tested by the user on return to the calling program. 

ICODE = 0: Normal return. 


f(xn) - f(xn-l) 
^n " ^n-1 
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APPENDIX D - Concluded 


Accuracy : The iteration process is continued until either of two convergence criteria 
are satisfied. These criteria are given as follows: 

If 

|f(xn) I = ei 

then 


f(xn) - Xn 
f(xn) 



( 1 ) 


and if 

|f(xn)| < ei 

then 




( 2 ) 


Reference : (a) Scarborough, James B.: Numerical Mathematical Analysis, Fourth ed., 
Johns Hopkins Press, 1958, p. 192. 

Storage : 137g. locations. 

Subroutine date: August 1, 1968. 
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APPENDIX E 


LANGLEY UBRARY SUBROUTINE FTLUP 

Language ; FORTRAN 

Purpose : Computes y = F(x) from a table of values using first or second order 
interpolation. An option to give y a constant value for any x is also provided. 

Use : CALL FTLUP (X, Y, M, N, VARI, YARD) 

X The name of the independent variable x 

Y The name of the dependent variable y = f(x) 

M The order of interpolation (an integer) 

I M = 0 for y a constant as explained in the NOTE below. 

M = 1 or 2. First or second order if VARI is strictly increasing 
(not equal). 

M = -1 or -2. First or second order if VARI is strictly decreasing 
(not equal). 

N The number of points in the table (an integer). 

VARI The name of a one -dimensional array which contains the N values of 

the independent variable. 

YARD The name of a one -dimensional array which contains the N values of 

the dependent variable. 

Note: VARD(I) corresponds to VARI(I) for I = 1, 2, . . . , N. For M = 0 or 

N <1, y = F(VARI(1)) for any value of x. The program extrapolates. 

Restrictions : All the numbers must be floating point. The values of the independent 
variable x in the table must be strictly increasing or strictly decreasing. The 
following arrays must be dimensioned by the calling program as indicated: VARI(N), 
VARD(N). 

Accuracy : A function of the order of interpolation used. 
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APPENDIX E - Concluded 


References : (a) Nielson, Kaj L.: Methods in Numerical Analysis. The Macmillan Co., 
C.1956, pp. 87-91. 

(b) Milne, William Edmund: Numerical Calculus. Princeton Univ. Press, 
C.1949, pp. 69-73. 

Storage : 430g locations. 

Error condition : If the VARI values are not in order, the subroutine will print "TABLE 
BELOW OUT OF ORDER FOR FTLUP AT POSITION xxx TABLE IS STORED IN 
LOCATION xxxxxx" (absolute). It then prints the contents of VARI and YARD and stops 
the program. 


Subroutine date: September 12, 1969. 



APPENDIX F 


LANGLEY LIBRARY SUBROUTINE DISCOT 


Language : FORTRAN 

Purpose: DISCOT performs single or double interpolation for continuous or discontinuous functions. 
Given a table of some function y with two independent variables, x and z, this subroutine performs 
Kxth- and K 2 th-order interpolation to calculate the dependent variable. In this subroutine all single- 
line functions are read in as two separate arrays and all multiline functions are read in as three 
separate arrays; that is, 


Xi 

(i= 1, 2, . 

. L) 

yj 

<N 

II 

. M) 


II 

. .,N) 


Use : CALL DISCOT (XA, ZA, TABX, TABY, TABZ, NC, NY, NZ, ANS) 


XA 

The X argument 


ZA 

The z argument (may be the same name as 

X on single lines) 

TABX 

A one-dimensional array of x values 


TABY 

A one-dimensional array of y values 


TABZ 

A one-dimensional array of z values 


NC 

A control word that consists of a sign (+ or -) 

and three digits. The control word is formed 


as follows: 


(1) If NX = NY, the sign is negative. If NX * NY, then NX is computed by DISCOT as 
NX = NY/Nz, and the sign is positive and may be omitted if desired. 

(2) A one in the hundreds position of the word indicates that no extrapolation occurs above 
Zmax' With a zero in this position, extrapolation occurs when z > z^ax- The zero 
may be omitted if desired. 

(3) A digit (1 to 7) in the tens position of the word indicates the order of interpolation in 
the x-direction. 

(4) A digit (1 to 7) in the units position of the word indicates the order of interpolation in 
the z-direction. 

NY The number of points in y array 

NZ The number of points in z array 

ANS The dependent variable y 



APPENDIX F - Continued 


The following programs will illustrate various ways to use DISCOT: 


CASE I: Given y = f(x) 

NY = 50 

NX (number of points in x array) = NY 
Extrapolation when z > Zjj^^x 
Second-order interpolation in x-direction 
No interpolation in z-direction 
Control word = -020 
DIMENSION TABX (50), TABY (50) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY 
READ (5,1) XA 

CALL DISCOT (XA, XA, TABX, TABY, TABY, -020, 50, 0, ANS) 

CASE II: Given y = f(x,z) 

NY = 800 
NZ = 10 

' NX = NY/NZ (computed by DISCOT) 

Extrapolation when z > Zjjj^x 
Linear interpolation in x-direction 
Linear interpolation in z-direction 
Control word =11 

DIMENSION TABX (800), TABY (800), TABZ (10) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY, TABZ 
•READ (5,1) XA, ZA 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, 11, 800, 10, ANS) 

CASE in: Given y = f(x,z) 

NY = 800 
NZ = 10 
NX = NY 

Extrapolation when z > 

Seventh-order interpolation in x-direction 
Third-order interpolation in z-direction 
Control word = -73 

DIMENSION TABX (800), TABY (800), TABZ (10) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY, TABZ 
READ (5,1) XA, ZA 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, -73, 800, 10, ANS) 

CASE W: Same as Case in with no extrapolation above Control word = -173 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, -173, 800, 10, ANS) 
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APPENDIX F - Continued 


Restrictions : See rule (5c) of section "Method” for restrictions on tabulating arrays and discontinuous 
functions. The order of interpolation in the x- and z-directions may be from 1 to 7. The following 
subprograms are used by DISCOT; UNS, DISSER, LAGRAN. 


Method : Lagrange's interpolation formula is used in both the x- and z-directions for interpolation. This 
method is explained in detail in reference (a) of this subroutine. For a search in either the x- or 
z-direction, the following rules are observed: 

(1) If X < Xj^, the routine chooses the following points for extrapolation: 

XpXg, . . .,Xj^^^ and y^,y 2 , . . ., yj^^^ 

(2) If X > Xjj, the routine chooses the following points for extrapolation: 

^n-k’*n-k+l’ • • -’^n yn-k’^n-k+l’ • • ^n 

(3) If X g Xf(, the routine chooses the following points for interpolation: 

When k is odd, 


X 

i 



When 

X 


k is even. 


. k’’'. k 1 ’ 
‘-2 




and 






y. k-y. k 

i -_ 1--+1 

2 2 



(4) If any of the subscripts in rule (3) become negative or greater than n (number of 
points), rules (1) and (2) apply. When discontinuous functions are tabulated, the indepen- 
dent variable at the point of discontinuity is repeated. 

(5) The subroutine will automatically examine the points selected before interpolation and if 
there is a discontinuity, the following rules apply. Let x^ and x^^j be the point of 
discontinuity. 

(a) If X £ x^, points previously chosen are modified for interpolation as shown: 

^d-k’^d-k+1’ • • -'^d yd-k’yd-k+l> ’ ' •> yd 

(b) If X > x^, points previously chosen are modified for interpolation as shown: 

^d+l’^d+2’ • • •^’'d+k yd+l’yd+2’ • • •> yd+k 

(c) When tabulating discontinuous functions, there must always be k + i points above 
and below the discontinuity in order to get proper interpolation. 

(6) When tabulating arrays for this subroutine, both independent variables must be in 
ascending order. 
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APPENDIX F - Concluded 


(7) In some engineering programs with many tables, it is quite desirable to read in one 
array of x values that could be used for all lines of a multiline function or different 
functions. Even though this situation is not always applicable, the subroutine has been 
written to handle it. This procedure not only saves much time in preparing tabular 
data, but also can save many locations previously used when every y coordinate had 
to have a corresponding x coordinate. Another additional feature that may be useful 
is the possibility of a multiline function with no extrapolation above the top line. 

Accuracy : A function of the order of interpolation used. 

Reference : (a) Nielsen, Kaj L.: Methods in Numerical Analysis. The Macmillan Co., c.1956. 

Storage : 555g locations. 

Subprograms used : UNS 40g locations. 

DISSER llOg locations. 

LAGRAN 55g locations. 

Subroutine date : August 1, 1968. 
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